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ABSTRACT 


Geolocation technology with the ability to locate an unknown beacon signal in three- 
dimensional space has been engrafted into numerous modern electronic systems. Indeed, 
the marketplace is anxious for more accurate and more accessible geolocation data. A 
primary limiting factor of the growth of geolocation systems is the stringent physical 


resource requirements needed for existing geolocation algorithms. 


Popular geolocation algorithms measure the time-of-arrival, time-difference-of- 
arrival, and frequency-difference-of-arrival of an incoming beacon signal from an 
unknown emitter at a given time. For these techniques, accurate solutions require a 
minimum of three airborne sensors; if available, a fourth sensor often significantly 
improves the accuracy. This resource requirement is excessive; we aim to relax it to two 
airborne sensors by applying a synthetic aperture technique. By fusing together data from 
multiple subsequent time samples, one can boost the overall resolution of the geolocation 


estimate. 


We propose using a series of geolocation measurements collected between two 
sensors according to a synthetic aperture model. System performance dependence on 
sensor velocity and aperture size is assessed. Additionally, a brief treatment of noise 
tolerance and estimation theory is given. Lastly, the overall feasibility of a synthetic- 


aperture-based geolocation algorithm is summarily addressed. 
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EXECUTIVE SUMMARY 


In recent years, the number of electronic systems using geolocation data has proliferated. 
A wide variety of growing applications is driving researchers to consider methods of 
extending coverage areas and boosting spatial resolution. Commonly used geolocation 
algorithms, including time-of-arrival (TOA), time-difference-of-arrival (TDOA), and 
frequency-difference-of-arrival (FDOA), require numerous airborne sensors—whether a 
constellation of satellites or a temporary network of aircraft—to deliver accurate results. 
In regions of operation with limited aerial view, such as canyons and valleys, it may be 
difficult to establish a line-of-sight communications link with many airborne sensors. The 
main research goal of this work is to reduce the number of required airborne sensors for 
accurate geolocation to two sensors. We seek to accomplish this using a synthetic 
aperture technique to combine subsequent geolocation measurements in time and develop 


a refined geolocation estimate. 


To lay an initial theoretical framework, we first consider the two- and three- 
dimensional solutions for the three primary algorithms under consideration: TOA, 
TDOA, and FDOA. The associated solution curves or surfaces are presented and 
discussed. To model these algorithms as random processes, we also introduce and define 
time- and frequency-domain noise models for various cases based on the sensor altitudes. 
For satellite-based sensors, we consider the effects of signal propagation through the 
ionosphere and the associated time delay and Doppler shift variations. For lower-altitude 
sensors onboard an unmanned aerial vehicle (UAV), we consider the time delay 
characteristics of lower-atmospheric propagation. For each scenario, an additive noise 


model is defined for TOA, TDOA, and FDOA measurements. 


The main thrust of our work is directed toward the implementation of a synthetic 
aperture algorithm with a two-sensor system. Accordingly, we first propose a 
homogeneous synthetic aperture TOA algorithm. It is considered homogeneous since 
both geolocation measurements collected at a given time step in the synthetic aperture are 
TOA measurements; there is only one type of measurement included in the algorithm. 


With TOA, we may collect one measurement per sensor in the network. In this case, we 


XV 


collect two TOA measurements per time step. Using a synthetic aperture of length JN , 


we systematically combine pairs of TOA measurements from subsequent time steps f, 


and ¢,,,. When the full duration of the synthetic aperture has passed, a refined estimate of 


the unknown emitter position is generated. Our numerical results for the synthetic 
aperture TOA algorithm show that, for accurate usage in a satellite-sensor system, future 
implementations will need to incorporate some noise mitigation techniques to improve 
the geolocation accuracy to useful levels. For UAV-based systems at altitude 10,000 feet, 
the results are more promising. Especially in low-noise scenarios, root-mean-square 
(RMS) estimation errors as low as 12.0 meters are achieved using the synthetic aperture 
technique. On the whole, we find that estimation error decreases as the number of 
synthetic aperture samples N increases; however, there is a diminishing return for large 
N as the RMS error asymptotically approaches zero meters. In such UAV-based 
systems, the synthetic aperture TOA algorithm is very promising. 


The second synthetic aperture algorithm we propose is a_ heterogeneous 
combination of one TDOA and one FDOA measurement at each time step in the aperture. 
For both the TDOA and FDOA geolocation techniques, we may collect one measurement 
per two sensors; thus, with two sensors, we have one TDOA and one FDOA 
measurement available at a given point in time. We again apply the synthetic aperture 
concept to synthesize data collected at subsequent steps in time and generate a refined 
geolocation estimate. A result of the complexity of the FDOA equation, we do not 
attempt an analytical solution to the TDOA-FDOA fusion problem and instead consider 
an iterative numerical solution. Our results for various conditions show that our proposed 


algorithm was limited by the excessive computational demands of our solution strategy. 
An iterative search through three-dimensional space scales as O(K \, where K is the 
number of discrete points in each dimension and o( -) represents traditional big O 


notation. We are limited in our ability to adequately sample the region of interest. 
Accordingly, the Nyquist sampling criterion bounds our ability to discern high-resolution 


spatial information. 
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Future implementations of a heterogeneous synthetic aperture TDOA-FDOA 
fusion algorithm will require reducing computational demand and addressing the Nyquist 
sampling limitation if it is to be suitable for high-accuracy applications. Some potential 
improvements to this end are suggested and discussed. Summarily, we find that a 
synthetic aperture TOA algorithm is both feasible and useful for reducing the physical 
resource requirement to two sensors for low-altitude UAV-based systems. The more 
complicated synthetic aperture TDOA-FDOA fusion algorithm will require a more 


optimized solution strategy to yield high-accuracy results. 
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I. INTRODUCTION 


A. BACKGROUND MOTIVATION 


In recent decades, terrestrial geolocation systems have shifted from the domain of 
academic research and niche military applications to that of the average household 
consumer. From the automotive industry to cellular phones, geolocation hardware and 
algorithms are now embedded onboard many everyday devices. Public access to the 
Global Positioning System (GPS) infrastructure has largely engrained into the popular 
mind that GPS is the only prominent geolocation system; however, the general topic of 
geolocation encompasses a much larger variety of systems, algorithms, and applications 
than GPS alone. Early competitors to GPS, including LORAN C, GLONASS, OMEGA, 
and GEOSTAR, all operated based on different geolocation methods and thus reported 


various performance characteristics [1]. 


Not all geolocation system applications are conventional or mundane. Researchers 
in Australia track the locations and frequency of lightning strikes in remote regions to 
improve wildfire containment strategies [2]. Geolocation data has been used extensively 
in seismologic research near the San Andreas Fault in California [3]. The E-911 system 
implemented in 2001 incorporates available GPS data with any proximate cellular 
infrastructure to offer emergency response services accurate locations of victims in 
distress [4]. In short, users from the military to the commercial marketplace are investing 
time and resources into developing geolocation-dependent systems. To support this 
growing demand, geolocation technologies must continue to progress and innovate. 
Without a readily available and robust technological infrastructure, geolocation 


applications may begin to stagnate. 


B. GEOLOCATION FUNDAMENTALS 


In general, a passive geolocation system consists of three conceptual components: 
the physical infrastructure, the geolocation measurement method, and the post-processing 
algorithm. The common physical scenario considered in the literature is that of a network 
of N airborne sensors with known position and velocity attempting to geolocate a single 
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emitter of unknown position and velocity. The emitter generates an electromagnetic 
beacon signal that propagates by line-of-sight (LOS) transmission to all sensors in view. 
Depending on the measurement method implemented by the system designers, the emitter 
and sensors may or may not be synchronized in time. The sensors may be satellites at 
various orbital altitudes or, as is the more recent trend in small-scale portable systems, the 


sensors may be unmanned aerial vehicles (UAVs). 


A variety of geolocation measurement methods are available, depending on 
system constraints. Angle-of-arrival (AOA) techniques use spatially-sensitive antenna 
arrays onboard each sensor to determine the direction from which an incoming signal 
arrives. With sufficient angular measurements, the sensors can determine the emitter 
location. The AOA method has largely been set aside in favor of others due to the 
complex array design necessary to achieve high angular resolution [5]. Time-of-arrival 
(TOA) techniques compute the absolute propagation time of the beacon signal from the 
emitter to each distinct sensor. From the propagation time, the distance from the emitter 
to each sensor is known; the emitter position can then be solved geometrically [6]. The 
GPS system employs the TOA measurement method [5]. Time-difference-of-arrival 
(TDOA) techniques compute the temporal difference between the beacon signal’s arrival 
time at multiple sensors. Here, no information is required of the emitter’s timing 
information; only the difference in arrival time between two sensors is of concern. 
Solving the TDOA measurement equation in three dimensions is more computationally 
intensive than AOA or TOA, but less a priori information is required of the emitter in 
order to geolocate [7]. Perhaps the most complicated method presented in this overview, 
the frequency-difference-of-arrival (FDOA) technique computes the emitter position 
based on the Doppler shifts seen at each individual sensor. Because of its Doppler 
dependency, the FDOA method requires some relative radial movement between the 
emitter and at least one sensor [8]. To be sure, other geolocation measurement methods 
exist, but the four discussed above cover a sufficiently wide topical breadth to suffice for 


introductory purposes. 


Much of the accuracy of the geolocation system can be attributed to the post- 
processing algorithms run following the measurement collection. There are typically two 
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stages of processing. First, an estimator algorithm fuses all the noisy measurements 
together and determines the best estimate of emitter position. For a noise-free 
environment, all measurements would coincide at a single spatial coordinate and 
estimation is unnecessary, as an exact solution is available. Unfortunately, all real 
implementations of geolocation systems suffer signal corruption due to noise. Resolving 
the measurement inconsistencies is traditionally done using one of two classes of 
estimators: least-squares (LS) estimators and maximum-likelihood (ML) estimators [7]. 
Each has its advantages and disadvantages, depending on the system specifics. Second, a 
tracking algorithm constructs predictions of future emitter position and velocity based on 
past measurement and estimation data. Many tracking filters have been presented in the 
literature, but two seem to have been established as field standards. The Extended 
Kalman Filter (EKF) attempts to linearize nonlinear measurements using Taylor 
expansions. For mild circumstances, the EKF works well in tracking applications; 
however, if the pertinent geolocation equations become severely nonlinear, EKF 
performance diminishes rapidly [9]. In response to some of the shortcomings of the EKF, 
the Unscented Kalman Filter (UKF) was developed as an alternative to avoid some of 
those linearization issues [9]. Other tracking filters have been proposed, to include 
various particle filters (see [10] and [11] for further discussion) and the Gaussian Mixture 
Measurement (GMM) filter [8]; however, one of the simpler Kalman Filter designs is 


often the more popular choice. 


C. RESEARCH GOAL 


Regardless of the specific geolocation method used, it is a general rule that N 
available sensors only allow geolocation in N -dimensional space with any degree of 
accuracy. So, in three dimensions, a constellation of at least three airborne sensors 
performing geolocation calculations is required to locate an unknown emitter. A system 
with only three sensors is said to be critically-determined. Many estimators allow for—or 
even prefer—more than three measurements in order to increase accuracy and spatial 
resolution; this is an over-determined condition. Alternatively, assumptions may be made 
about the physical environment; assuming a known emitter altitude, for instance, narrows 


the three-dimensional search to two dimensions, namely the horizontal plane at the 
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known altitude. This, in turn, relaxes the requirement on number of measurements and 


sensors [12]. 


For many conventional operating scenarios, the minimum resource requirement of 
three visible airborne sensors is a nonissue. Often either geographic landscapes offer a 
sufficiently wide aerial field-of-view (FOV) to accommodate LOS transmissions with 
three or more sensors, or the end user is willing to suffer gaps in geolocation coverage 
areas. Accordingly, little effort to date has been directed toward assessing techniques to 
relax the resource requirement; however, it is not difficult to imagine a situation where 
the three-sensor requirement is unacceptable. For instance, one can easily imagine an 
urgent search and rescue application that demands accurate geolocation fixes in 
mountainous terrain with limited FOV and limited sensor availability. Currently, when a 
system is below this critically-determined threshold, it will fail to provide geolocation 
information with any helpful degree of accuracy; until the resource requirement can be 
met, the emitter must bear the burden of estimating its own position. Depending on the 
application, local onboard solutions may be available. For instance, automobile 
navigation systems couple inertial sensor data and accurate road maps to help interpolate 
an automobile’s position between intermittent GPS fixes. The inertial sensors provide 
information on course changes, while the road maps restrict possible vehicle movements 
along known paths. Together with recent GPS data, this additional local information 
helps offer a seamless position estimate in real time [13]. Still, if size and power 
requirements for the emitter are strict, inertial sensor suites and significant onboard 
processing power is an excessive demand. Rather than push the requirement complexity 
from the sensors to the end user’s emitter, a new accurate geolocation technique with 


only two required airborne sensors is preferable. 


In 2007, Fletcher, Ristic, and Musicki proposed using the TDOA measurements 
from only two moving UAVs to geolocate an emitter [14]. They proposed a recursive 
estimator, which updates the estimated emitter position following each successive 
measurement rather than entirely recalculating the emitter position based on previous 
measurement history. They examined both the EKF and UKF as potential tracking 
algorithms for both a stationary and moving emitter. Each filter performed accurately and 
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converged on the Cramér-Rao Lower Bound (CRLB)—the theoretical lower limit for 
mean-squared error (MSE) of any unbiased estimator [14]. Their result introduced a new 
concept in reducing the physical resource requirement for accurate geolocation. One can 
substitute temporal resources for physical resources. An under-determined network of 
sensors might be able to geolocate an emitter by fusing together a series of successive 
measurements in time. Moreover, each additional measurement may be incorporated to 
improve the spatial resolution of the algorithm. This functional concept has been 
implemented for years by synthetic aperture radars (SAR), which combine many 
successive returns from a low-resolution radar to generate a very high-resolution image. 
As with SAR, each successive measurement increases the accuracy of the estimated 


emitter position. 


It should be noted that Fletcher, Ristic, and Musicki collected one TDOA 
measurement per time step [14]. This homogeneous measurement technique may have 
advantages in computational efficiency, but it is possible to extract more measurement 
data from the two available sensors. In this research, we propose and investigate a two- 
sensor geolocation algorithm that collects at each time step one TDOA measurement and 
one FDOA measurement. This heterogeneous technique fuses together the TDOA and 
FDOA measurement types. Incorporating FDOA computations undoubtedly increases the 
complexity of the solution approach and raises the demand for processing power; 
nevertheless, doubling the number of measurements per time period may offer 
worthwhile improvement in overall accuracy. As a preliminary step toward pursuing the 
TDOA/FDOA data fusion, we first consider the simpler case of fusing together two TOA 


measurements as a proof of concept. 
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I. TIME-OF-ARRIVAL 


A. GOVERNING EQUATION 


On the heels of the widespread success of GPS, TOA-based geolocation 
algorithms have risen in prominence. The solution method is simple when compared with 
methods such as TDOA and FDOA. Still, achieving simplicity in an algorithm often 
demands more numerous or detailed algorithm inputs. One key aspect of the TOA 
algorithm is that absolute timing information is necessary. Namely, when each individual 


sensor i receives the emitter signal, it must be able to determine the absolute one-way 


signal propagation time. This is the TOA measurement 7,;. Multiplying by the 


propagation speed Cy), we relate propagation time to the range between the emitter and 


sensor i. Assuming a noise-free environment and a constellation of N TOA sensors, the 


emitter position must fall within the set of all points that satisfy 
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=ct,, i=1,2,...N, (1) 


where Fr, and F; are the range vectors from the coordinate system origin to the emitter 


L 


and to sensor i, respectively. The vector 7, is simply the difference between vectors F, 


and 7. Unless specified otherwise, bold-faced font hereafter indicates a vector. 
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Additionally, the } : 1 operator hereafter denotes the Euclidean norm of a vector. 


B. TWO-DIMENSIONAL SOLUTION 
Consider a two-dimensional scenario as depicted in Figure 1. The two sensors 
labeled 5, and S, attempt to geolocate an emitter of unknown position labeled e . In this 


TOA geolocation scenario, there are no restrictions on sensor or emitter movement. The 
specific choices of sensor position, emitter position, and units are arbitrary, so convenient 


values are chosen for presentation’s sake. The sensors are located at the known Cartesian 


coordinates (—5,0) km and (5,0) km; the unknown emitter is located at (1,7) km. 


y-axis (km) 
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x-axis (km) 
Figure 1. Two TOA sensors 5, and S, attempt to locate an emitter of unknown 
position. 


As in the general case for TOA geolocation, the emitter generates an 
omnidirectional beacon signal that propagates toward each sensor. Upon signal arrival, 


each sensor determines the total transmission time of the signal from emitter to that 


respective sensor. The TOA measurement collected by sensor 5S, in turn gives the radial 
range C,Z, to the emitter position. This restricts the emitter position to the set of all points 


on the circle with radius C,7,and centered at the corresponding sensor. This circle may be 
referred to as a constant-TOA curve; an emitter located at any point on the curve 


generates the same TOA measurement value at sensor 5; . 


For the two-sensor scenario ( N = 2 ), the resultant circular emitter position curves 
are shown in Figure 2. Encouragingly, both circles intersect at the emitter location; 


however, there is a point of ambiguity where both circles intersect at a point other than 
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the emitter location. As such, even in two-dimensional space, two TOA measurements 
alone are insufficient to unambiguously locate the emitter. More measurements are 


required to give a singular solution. 
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Figure 2. Two TOA emitter position curves are color-coded to match the corresponding 
sensor. 


To solve the ambiguity problem, an additional sensor 5, is positioned at (—5,10) 


km and calculates an additional TOA measurement based on the received emitter signal. 


As with the other sensors, we solve the TOA measurement for the emitter position curve. 


The additional position circle from sensor 5; is shown in Figure 3 atop the previous 


measurements from sensors 5S, and S,. 
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Figure 3. Three TOA emitter position curves unambiguously geolocate the unknown 
emitter. 


Here, the third sensor has resolved the ambiguity issue. There are still multiple 
points where two of the three position curves intersect, but there is only one single 
location where each position curve intersects all others. So, with three sensors we have 
accurately and unambiguously located the unknown emitter in a two-dimensional, 


noiseless environment using circular TOA measurements. 


Still, even without noise corruption, there are situations in which no number of 
TOA sensors can unambiguously located an emitter. Consider the scenario wherein 


sensor 5; is placed at (0, 0) km and calculates a TOA measurement as expected. The 


resulting emitter position curves are plotted in Figure 4. 
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Figure 4. | Three TOA emitter position curves unsuccessfully attempt to resolve emitter 
location ambiguity. 


Clearly, the addition of the third sensor does not resolve points of ambiguity 
under all circumstances. In fact, for any collinear arrangement of TOA sensors, we are 
always left with an ambiguous solution, irrespective of the number of sensors present. 
This is a point of great concern; the physical arrangement of geolocation sensors can have 
a severe impact on the performance and capability of the geolocation system. Moving 
forward, it is important to characterize the capabilities of various geolocation techniques 


with respect to physical orientation and movement of the sensor constellation. 


C. THREE-DIMENSIONAL SOLUTION 


Extension of the TOA solution to three dimensions is not complicated. Rather 


than a two-dimensional circle, the solution to Equation 1 is now a sphere centered around 


the i-th sensor position with radius C,7,. The emitter’s unknown position is necessaril 
p 0°71 Pp y 
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restricted to the set of points lying on the resultant spherical TOA surface. A second TOA 
measurement generates another spherical emitter surface. In general, the intersection of 
any two spheres is either a single point or a circle. If the spheres intersect tangentially, 
there will be only a single point of intersection, but this condition is effectively a 
theoretical singularity that does not occur in practice. Alternatively, two spheres may 


intersect to form a circular intersection curve. This principle is illustrated in Figure 5. 
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Figure 5. | Two TOA emitter position surfaces restrict the potential emitter location to a 
circular curve. 


Here, two arbitrarily positioned TOA sensors are placed at the Cartesian 
coordinates (—5,0,0) km and (5,0,0) km and depicted as red markers. The emitter is 
located at (-4, —2,6) km and depicted as a black marker. The sensors receive the TOA 


beacon signal and calculate the associated spherical emitter surfaces. The circular 


intersection of the two spheres is shown as a dark blue curve. Here, since the TOA 
12 


sensors are both located directly on the x-axis, the circular intersection curve is contained 
in a vertical plane perpendicular to the x-axis. In fact, for the two-sensor condition, the 
circular intersection curve is always contained within a plane perpendicular to the sensor 
axis, regardless of sensor orientation. Notice further that combining two TOA sensor 
measurements effectively brings us back to a two-dimensional scenario as discussed 
previously. It follows, then, that an additional TOA sensor positioned away from the 
existing sensor axis further restricts the emitter position to two potential points. Finally, a 
fourth TOA sensor positioned away from all existing sensor axes solves for the 
unambiguous emitter position at a single three-dimensional coordinate. We add the 
caveat that additional sensors are placed away from existing sensor axes to avoid the 
insoluble ambiguity problem that arises from all sensors being placed on the same sensor 
axis, as discussed previously. Altogether, extension of the TOA solution method into a 
third dimension increases the sensor requirement by one sensor in order to preserve the 


same degree of accuracy. 
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HW. TIME-DIFFERENCE-OF-ARRIVAL 


A. GOVERNING EQUATION 


A less restrictive alternative to TOA is TDOA, which does not require knowledge 
of the absolute propagation time of the emitter’s beacon signal. Instead, the TDOA 
geolocation technique uses highly synchronized clocks onboard all devices in the sensor 
constellation to monitor the relative difference in arrival time between sensors. Each 


TDOA emitter position curve or surface is defined as the set of all points satisfying 


rT, ~ i, j =1,2,...,N, i# J, (2) 


e 




















r =r,|=¢t;;5 


where the time-difference-of-arrival measurement given by 
t,,=(t,-%)-(t,-%)=t,-7; (3) 
may be defined in reference to the absolute time of transmission 7, or simply as the 
relative difference in arrival time at sensors s, and s, , denoted t,—7,. Note that the 
TDOA measurement 7,, may take on any real value, either positive or negative. From the 


study of analytical conics, recall that the set of all points such that the difference in 
distance to two foci is constant is a hyperbolic curve [15]. Accordingly, if the problem is 
cast in two-dimensional space, the particular solution curve for Equation 2 is a hyperbola; 
if it is cast in three-dimensional space, the particular solution surface for Equation 2 is a 


hyperboloid of two sheets. Sensors s; and s, are located along the major axis of the 


hyperbola or hyperboloid and function as the hyperbolic foci. 


B. TWO-DIMENSIONAL SOLUTION 


In the two-dimensional x-y plane, two sensors of known position attempt to locate 


the position of an unknown emitter using TDOA measurements. In the following section, 
sensors S, and S, are located at the generic Cartesian coordinates (% y,) and (%5 Yo) ; 
respectively. It is again assumed that the environment is noise-free and all measurement 


values are exact. 
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In this scenario, the emitter begins transmitting an omnidirectional beacon signal 


at some unknown time 7). This signal arrives at sensor S, at time 7, and at sensor S, at 


time T,. The lack of knowledge of the originating time 7, prohibits the use of TOA 


measurements in this case; however, a TDOA measurement may still be collected. 
Comparing the difference in signal arrival time at each sensor, we solve (3) for the single 


TDOA measurement 7,,=7,—7,. With a known TDOA measurement value, we may 


solve (2) analytically. 


The following derivation is an algebraic solution to (2) based on the proof given 


in [16], with some appropriate modifications. We solve for the set of all Cartesian 


coordinates (og y) that satisfy the TDOA measurement. Using the distance formula, we 


rewrite the TDOA equation 














r -r,|- r, —A |= Cota, (4) 


e 


as 


W(x-%) +(y-y) —(x-x,) +(y-y,) = CT: (5) 


We next orient our coordinate systems such that sensors S, and S, are located at the 


Cartesian coordinates (—d,0) and (d,0), respectively. The problem then reduces to 


i(x-dyY +y? -f(x+d) +y’ = C57 91° (6) 


Rearranging terms and squaring both sides of the equation gives 
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| (say +y") =(ora 4 (xtra) +5"). (7) 


We subsequently expand to get 





x -2dxt+d*+y? =(¢,t) + CT VX? +2dx+d?t+y +x +2dx+d>+y. (8) 
Collecting all linear terms and constants on the left-hand side of the equation, we are left 


with 
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—Adx—(c,t) = 2(cytm) Vx +2dx+d?+y’. (9) 


To condense the notation, we introduce the parameter 


Cota 
a=——., 10 
5 (10) 


which carries the same sign as the TDOA measurement 7,,. Substituting (10) into (9) 


gives 


—dx—a’ =a,Jx’+2dx+d’+y°. (11) 

We then square both sides of the equation and expand to get 
d°x? +2da°x+a* =a? (x? +2dx+d’+y’). (12) 

Note that the sign information associated with parameter @—and thus the TDOA 
measurement T,,—is now lost from the squaring operation. We shall revisit this issue 
below. Collecting terms and factoring leaves 

(d?-a’)x? -a’y’ =a’ (d’ -a’). (13) 
We then divide through by a’ (d > a’) to give 


2 y? 
a are =1. (14) 


Introduction of the traditional hyperbolic parameter 


b=Va@?-@ (15) 


into (14) gives the recognizable form of a hyperbola with foci lying along the x-axis, 





namely 


x y’ 

—--—=1. 16 

og (16) 
To recast the hyperbola in terms of quantifiable TDOA parameters, we substitute out the 


intermediate parameters a and b, given by (10) and (15), respectively, to get 


Age = Ay? 
(cf) 4d? —(€97;) 
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=I; (17) 


Note that (17) is fully expressed in terms of the spatial Cartesian coordinates 


(~y), the known propagation speed C,, the known spatial parameter d, and the 


calculated TDOA measurement 1,,. Altogether, we have sufficient information to 
localize the unknown emitter to a two-dimensional hyperbolic curve. 

One final point requires discussion. The hyperbolic equation in (17) is satisfied by 
two branches. Let us term these the positive and negative branches. The positive branch 
is that half of the hyperbola nearest in space to sensor s,, and the negative branch is that 
half nearest in space to sensor s,. Since spatial distance directly correlates to time delay, 
we know that an emitter located at any point on the positive branch corresponds to a 
TDOA measurement 7,, =7,—7, > 0. Conversely, an emitter located at any point on the 


negative branch must correspond to 7,,<0. We noted above that the sign information 


associated with the TDOA measurement 7,, was lost in (12) due to the squaring 


operation. We may now retroactively go back and rule out one of the two hyperbolic 
branches. In short, choose the single branch that corresponds to the sign of the received 


TDOA measurement and discard the other branch as an invalid solution. 

To illustrate the TDOA solution curve, we consider an example two-sensor 
TDOA problem depicted in Figure 6. The sensor and emitter positions in Figure 6 are 
arbitrarily chosen to aid visualization; the sensors 5, and S, are positioned at coordinates 
(—5,0) km and (5,0) km, and the emitter is positioned at the coordinates (4,7) km. The 
sensors are unaware of the emitter location prior to the arrival of the TDOA signal. 

We assume each sensor has preexisting knowledge of the spatial coordinates 


of all other sensors in the constellation. From the problem statement, the 


spatial parameter d=5 km is known. The speed of propagation of the emitter 
signal is taken to be the free-space value 2.998x10° m/s. Assuming ideal LOS 


signal transmission, the TDOA measurement between sensors S, and S, upon arrival is 


r 


e 


























) 6 =-14.44us. With only knowledge of these three parameters and 


Ty =( el 
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measurements, the sensors can solve for a particular solution to Equation 17. Figure 7 is a 
two-dimensional plot of the corresponding hyperbolic emitter curve. The solid blue trace 
denotes the valid solution branch; the dashed blue trace is invalid according to a sign 


difference with the TDOA measurement and may be discarded. 
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Figure 6. Two TDOA sensors attempt to locate an emitter of unknown position. 
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Figure 7. A single TDOA solution curve localizes the emitter to a hyperbolic curve. 
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As with the TOA technique, inclusion of a second measurement helps narrow the 


search for the emitter. Building on the previous example, we consider a third TDOA 
sensor 5; positioned at coordinates (-3,10) km. With three sensors in the constellation, 
there are three potential sensor pairs. As a first step, we include the measurement 
generated by the arrival time difference between sensors S, and S,. Since this sensor pair 


is no longer solely located along the x-axis or centered at the origin, the hyperbolic 
solution derived above must be rotated and translated to accommodate the new sensor 


orientation. 


To compute the required rotation matrix, we begin with the dot product 


expression 
r,,-% = r,,|]¥||cos(A.,, ), (18) 
where 7;, is the position vector originating at sensor s, and terminating at sensor s;, X 


is the positive x-direction unit vector, and @,,, is the angle of rotation between these two 


rot 

















vectors. Understanding that |x||=1, (18) may be rearranged as 
cos(6,,) ==. (19) 
[rs 


We temporarily restrict the rotation angle to the range (0°,180°) and solve: 


6, =arccos 7) (20) 


Irs 
To include rotation angles in the range (—180°,180°) , we must account for the inherent 


quadrant ambiguity of the arccosine function. Consider the signum function that is 


piecewise-defined as 


-l, a<0O 
sgn(a)=) 0, a@=0. (21) 
1, a>OoO 
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Note that the expression sgn(r;, : y) is positive for all points in the first and second 


quadrants and negative for all points in the third and fourth quadrants. We may 
accordingly compensate for the aforementioned arccosine quadrant ambiguity by writing 


31 
In| 


6, =sgn(r;, Saco 2 =| (22) 


Based on the rotation angle 0.,,, we define the rotation matrix [17] 


rot ? 


ees —sin(6,,, : (23) 


sin (6..,, ) COS (O,ur ) 


We further define the translation matrix 


> 


r, +7, - 
2 

r+n \ 

[ 2 }s 


which specifies the coordinates of the midpoint between sensors S, ands, in column 


T= (24) 


2 


vector format. The rotated and translated hyperbola coordinates fx, y’) are then given by 


x x 
=| Jer (25) 
y y 


Following the prescribed rotation and translation procedure, the second TDOA 
emitter curve generated from the measurement between sensors 5, and S$, is shown in 


Figure 8. This second TDOA hyperbola is depicted as a red trace; the original TDOA 


hyperbola from Figure 7 is reproduced as a blue trace, sans the invalid branch. For further 
assurance of an accurate solution, one might use the third sensor pair S, and S,. In 


Figure 9, the third TDOA solution curve is overlaid on the previous two TDOA curves. 
For visualization, the third solution curve is depicted as a red trace and the previous 


curves are reproduced as blue traces. 
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Figure 8. Two TDOA solution curves help locate an unknown emitter. 
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Figure 9. Three TDOA solution curves help locate an unknown emitter. 
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If more sensors are available, the corresponding measurements may be 
incorporated into the solution to increase estimation accuracy. In general, the number of 
possible sensor pair combinations for an N -sensor constellation is 

N N! 
o-(3 ae ov 
As a general trend, one would expect geolocation ambiguity to decrease as the number of 


TDOA measurements increases; yet, this does not always hold. 


Unfortunately, the TDOA geolocation technique is subject to the same 
orientation-dependent restrictions as the TOA technique. If the TDOA sensors in a given 
constellation are arranged in a collinear orientation, the system will never be able to 
unambiguously resolve the emitter location. In two dimensions, all TDOA hyperbolas 
from collinear sensors intersect at both the true emitter location and its reflection point 


about the sensor axis. This inherent orientation-dependency is depicted in Figure 10. 
20 
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Figure 10. | Four TDOA sensors unsuccessfully attempt to resolve geolocation ambiguity. 
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Here, four sensors are non-uniformly spaced along the x-axis, and three of the six 


possible TDOA hyperbolas are plotted. All emitter curves intersect at both the emitter 


location and the point of ambiguity a, at coordinates (4, -7) km. As such, the network 
can never successfully resolve the emitter location. 


Intuitively, it makes sense that collinear sensors entail an element of inherent 
ambiguity. First, note that the TDOA (and TOA) solution curve is always symmetric with 
respect to the sensor axis. For the TDOA case, the sensor axis is simply the major axis of 
the hyperbola. Second, note that all sensors lie along the same linear sensor axis (the x- 
axis in Figure 10). Together, these observations leave the sensors without any ability to 
discern knowledge regarding directions orthogonal to the sensor axis. In two dimensions, 


only two unit vectors satisfy this orthogonality criterion (+y and —y in Figure 10); one 
unit vector corresponds to the true emitter position (+y) and one to the point of 
ambiguity (—y ). As we shall see, expanding the collinear ambiguity problem into three 


dimensions is crippling. In three dimensions, there is now an infinitude of unit vectors 
that are orthogonal to the sensor axis. As the TDOA problem extends into higher- 
dimensional scenarios, one must be acutely aware of the spatial dependencies and 


limitations. 


Cc; THREE-DIMENSIONAL SOLUTION 


The three-dimensional solution to the TDOA governing equation (see Equation 2) 
is, predictably, analogous to the two-dimensional solution. Whereas in two dimensions 


the emitter curve is a hyperbola, in three dimensions the emitter surface is a hyperboloid 


of two sheets. Consider a generic scenario with two TDOA sensors 5, and S, located at 
coordinates (—d,0,0) km and (d,0,0) km, respectively. Assume the parameter d >0. An 
emitter of unknown location generates a TDOA beacon signal that arrives at sensors 5S, 


and Ss, with time difference of 7,,. It can be shown algebraically (see Appendix for 
details) that the associated emitter surface is the set of all three-dimensional points 


(x y, z) that satisfy the equation 
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4x? 4y? 42° 
<--> _, -___4_,=1. (27) 
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As in two dimensions, we refine this solution by selecting only the valid sheet according 
to the sign information of the TDOA measurement. The invalid sheet—corresponding to 


a TDOA measurement of the opposite sign—is discarded. 


For illustration, we contrive an arbitrary scenario where d=5 km and the 


unknown emitter is located at (-6, 4,5) km. The associated TDOA measurement value 
is calculated by sensors s, and S, to be 7,, =20.84 us. A three-dimensional plot of the 
hyperboloidal TDOA emitter surface is presented in Figure 11. In the figure, sensors 5S, 


and S, are depicted as red markers; the emitter is depicted as a black marker. 
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Figure 11. | Two sensors generate one TDOA surface in three dimensions. 
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Including a second TDOA measurement further refines the emitter location. For 
an unambiguous TDOA solution, it is commonly noted that four or more sensors are 
required [12]. The use of three sensors at best refines the emitter location to several 
discrete locations. In this case, a priori information such as altitude or terrain information 
may be used to resolve the emitter position, but four or more sensors are necessary to 
ensure an unambiguous solution. Still, this conclusion bears the previously discussed 
caveat regarding collinear sensor orientation. In the collinear sensor case, the TDOA 
geolocation technique can only specify a circular curve on which the emitter must lie. 
Again, an analogy with the two-dimension ambiguity case is helpful. Referring to Figure 
10, this circular ambiguity curve is simply the orbit traced by the points e and a when 
revolved around the sensor axis. It is therefore imperative that a TDOA sensor 
constellation avoid scenarios in which multiple sensors are arranged in collinear fashion. 
The associated TDOA measurements will unavoidably be redundant and the solutions 
ambiguous. Again, we note that the spatial orientation of the sensor constellation acutely 


affects the accuracy and functionality of the TDOA geolocation technique. 
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IV. FREQUENCY-DIFFERENCE-OF-ARRIVAL 


A. DOPPLER EFFECT 


The FDOA geolocation algorithm operates on the basic principle of the Doppler 
Effect. Consider a simple transmitter and receiver model as illustrated in the Figure 12 


system diagram. The transmitter generates a cosinusoidal signal with center frequency f, 
that propagates through free space at the speed of light c, and is received by the receiver 


(classically termed “the observer’). 


Transmitter, v, Receiver, v, 


Spectrum Analyzer 





Figure 12. _—_— A transmitter/receiver model experiences a Doppler shift. 


When both devices are stationary, the received frequency is the same as the 


transmitted frequency, namely f,. Yet if the one or both of the devices are in motion with 


respect to some common reference, the transmitter will have a nonzero relative velocity 
with respect to the receiver. The obvious exception, of course, is the case when the 
transmitter and receiver have identical absolute velocities; subsequently, the relative 
velocity is zero. Mathematically, we define the vectorized transmitter velocity relative to 


the receiver as 
Vip = V, —Vi> (28) 


where V, and ¥, are the respective transmitter and receiver velocities in vector form. 


When the relative transmitter velocity contains a nonzero radial component—that is, the 
transmitter is moving either toward or away from the receiver—the receiver observes a 
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shift in the observed frequency f,,.. If the transmitter is moving toward the receiver, the 
observed waveform appears compressed at the receiver such that f,,, > fy. Conversely, if 


the transmitter is moving away from the receiver, the observed waveform appears 


elongated at the receiver such that f,,.< f >. To quantify the frequency shift, we first 


define the unit vector 


2 nr, 
od 


orl’ 7 








where % and I, are vectors from the origin to the transmitter position and receiver 


position in space, respectively. Thus, the unit vector i, points from the receiver position 


to the transmitter position. To then extract the radial component of the relative transmitter 


velocity, we compute the scalar speed 


vi=y -i. (30) 
We define the associated Doppler shift as the difference between the original oscillator 


frequency and the observed frequency, namely 


A fp = Sots — So (31) 
We then relate the Doppler shift to the radial component of the relative transmitter 


velocity via the fundamental Doppler relationship [8] 


Af, =2(vi,). (32) 


Note again that, to achieve a nonzero Doppler shift, we require a nonzero radial 


component of the transmitter velocity relative to the receiver velocity. 


B. GOVERNING EQUATION 


The basic premise of the FDOA geolocation algorithm is to use the observed 
Doppler effect on an unknown emitter’s transmitted beacon signal to geolocate it in 
space. In the following section, we aim to derive a generalized FDOA governing equation 
based on the method proposed by Musicki and Koch [8]. Consider the generalized two- 


dimensional scenario depicted in Figure 13 to help illustrate the problem. For a single 
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FDOA measurement, two sensors are required with narrow frequency resolution. Sensors 
S, ands, move with the respective absolute velocities ¥, and v¥,. Assume the 
propagation medium is isotropic and ideal and the system is noiseless. The emitter e 


moves with absolute velocity V, and generates an omnidirectional FDOA beacon signal 


at the known oscillator center frequency f). The signal arrives at sensors Ss, and S, with 


observed center frequencies f, and f,. 


y-axis (km) 
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Figure 13. | A two-sensor FDOA geolocation scenario is depicted. 


From the observed frequencies measured at the two sensors, the FDOA 


measurement is calculated as 


Af =hi-Sfi- (33) 


Since f, and f, are both observed frequencies, (31) enables us to write 
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=A tho (34) 


and 
f, =A + ho, (35) 
where Af, and Af, are the Doppler shifts associated with the arriving signal at sensors 


S, and S,. Thus, the FDOA measurement may be expressed as the difference in Doppler 


shifts between the two sensors; substitute (34) and (35) into (33) to give 


Afn =(Af, + hy)-(AA + hy) =A NAL. (36) 
To translate this frequency information into position information, we must first 


express the FDOA measurement in terms of relative velocities. Of particular interest, two 


relative emitter velocities are defined with respect to sensors S, and S,: 
Vo =, Vis (37) 
Veo =V, —V». (38) 
According the generic form of (29), we next define the two similar unit vectors 


> r,—y, r 
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Referencing (30) and (32), we may now express (36) as 


Af, = Jf -vi,) 


Co 


_ to 


Co 


(arte ¥idg): (41) 
Equation 41 may be considered the fundamental governing equation for the FDOA 
algorithm. Here, the FDOA measurement between sensors 5S, and S, is linked with the 
position information contained in the unit vectors - and i... With all other system 
variables assumed to be known, any received FDOA measurement A/f,, can only result 


from a particular set of combinations of i, and bs This set of unit vector combinations 
30 


describes the constant-FDOA curve. Each pair of associated unit vectors in the set may 
be thought of as a pair of rays that identify a single point on the constant-FDOA curve. 
The two rays originate at the known sensor locations and extend in the direction specified 
by the associated unit vectors. We solve for the intersection of these two rays using a 


three-dimensional triangulation method. 


C. TRIANGULATION 


A brief treatment of the triangulation method used in this research will likely be 
helpful. As may be discerned by cursory inspection, the FDOA governing equation 
presented in (41) requires a nontrivial solution strategy. To provide helpful illustrations 
without excessive computational demands, we consider a parametric solution. Such a 
solution lends itself to discretized rendering in computational software environments 
(e.g, MATLAB"). Accordingly, we consider the problem of triangulation in parametric 


terms. 


The process of triangulation may be stated as follows. We attempt to determine 
the singular location in N-dimensional space at which two rays of known origination 
intersect, provided that they do indeed intersect. To begin the mathematical solution, a 


ray is defined in vector form as 

P(t)=it,+n, 1,20, (42) 
where 7, describes the ray origination point with respect to the coordinate system origin, 
i, is a unit vector parallel to the direction of the ray’s extension from the origination 


point, and f, is a nonnegative parameter. The ray p, is known to intersect the similarly- 


defined ray 

pp(t,)=it,+n, 1,20. (43) 
At the point of intersection, p, (1) =P, (1.} . Using primed notation to specify this point 
of interest, we have 


it, +7, =it, +n. (44) 
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Rearranging gives 


Uy 


r,—r, =i, -i,t,. (45) 
Let us define Ar,,=r,—r,, which may be decomposed into its three-dimensional 


Cartesian components by 





Ar,, = An, x+ Ar, ¥ + An, ,z. (46) 


We likewise decompose the unit vectors i, and i, as 


> 


~. 


hy Phy yt hae, (47) 


na 


tL, =b,Xth,Vth,z (48) 


Equation 46 may then be expanded into the system of linear equations 


An, = Ae ty i), (49) 
aa t& 

Any =thh, -t,b, (50) 
> aA 

Br Ty Big bs loys (51) 


To solve for the two unknown parameters 1, and 1,', only two of the three 
equations in the system are required. Rearranging (50) gives 
' 1 wn 
fo | (|, Ly - Ar}. (52) 
by 
Substituting (52) into (49) then gives 
LEAD | aa na 
An, =th, -| =— (" hy - An, |(é.): (53) 
2y 


Solving, we have 


Any. -{}: Jax, 
Uy l y 
ae (54) 


Thus, the two rays defined above intersect at the point specified by p, (1) , where ¢,' is 
computed as specified in (54). 


This parametric solution to the triangulation problem presumes the foreknowledge 
of both unit vectors that describe the two intersecting rays. In the FDOA problem, there is 
a pair of unit vectors describing each point along the constant-FDOA solution curve. 
Hence, to find all points on the continuous FDOA curve, one must perform this 
parametric triangulation an infinite number of times. This impossibility is obviously 


avoided. Rather, a sampled solution of the FDOA solution curve shall be computed. 


D. TWO-DIMENSIONAL SOLUTION 


The complexity of the FDOA solution curve is readily apparent even in the 
simplified two-dimensional case. We shall consider several scenarios in this section to 
illustrate the spatial dependencies of the FDOA algorithm. In particular, we demonstrate 


the impact of sensor motion relative to the sensor axis. Consider a generic scenario with 


two sensors 5; and S, moving with velocities Vv, and Vv, and located at coordinates 


(—d,0,0) and (d,0,0), respectively. We illustrate this setup for d =5 km in Figure 14. 
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Figure 14. | Two sensors move with respect to the sensor axis. 
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We shall consider two limiting conditions and one transitional condition. First, we 
shall consider the scenario where sensor velocities are identical and parallel to the sensor 
axis. Second, we shall consider the scenario where sensor velocities are identical and 
orthogonal to the sensor axis. Third, we shall consider a transitional scenario where 


sensor velocities are identical and oriented at 60° relative to the sensor axis. 


1. Parallel Motion Condition 


Consider the case where both sensors 5, and S, have identical velocities directed 


parallel to the sensor axis (here, the +x direction). The magnitude of the velocity is 


inconsequential; however, for illustration’s sake we choose the arbitrary sensor veloci 
> 


of +10% m/s. Sensors s, and s, are located at (—5,0) and (5,0) km, respectively. The 
emitter location—tecall this information is unknown to the system a priori—is chosen to 


be (6,4) km. The associated FDOA solution curve is presented in Figure 15. Note that 


the parallel motion condition yields a continuous constant-FDOA solution curve. 
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Figure 15. | Two FDOA sensors geolocate in the parallel motion condition. 
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We stated above that the magnitude of the sensor velocities was inconsequential. 
Indeed, if the FDOA solution curve is computed for a sensor velocity of +100 m/s, one 
finds an identical solution to the +10% m/s case. Certainly, with the tenfold increase in 
sensor velocity, the absolute value of the FDOA measurement increases proportionally; 
however, the FDOA curve is generated by triangulation of intersecting unit vectors. The 
issue of concern (see Equation 41) is not the absolute values of the sensor velocities but 
rather the projection of the sensor velocity vectors onto the associated unit vectors via the 
dot product. The FDOA solution curve is bound up with the directional information of 
the sensor velocities relative to the sensor axis. Thus, if the velocities are directed away 


from the sensor axis, we see a dramatic shift in form for the solution curve. 


2. Orthogonal Motion Condition 


To that end, we next consider the case where sensors 5, and S, have identical 
velocities directed orthogonally to the sensor axis, namely the +y direction. The 
magnitude of the velocity is again inconsequential, but we choose the arbitrary sensor 


velocity of +105 m/s. The emitter is at coordinates (6,4) km. See Figure 16 for the 


associated FDOA solution curve. Here, we see that the orthogonal motion condition 


yields a non-continuous, multiple-branch FDOA solution curve. 


The complexity of the FDOA solution become apparent. Whereas the TDOA 
algorithm yielded a hyperbola of varying parameter values, the FDOA algorithm may 
yield a disjointed, piecewise solution depending on the sensor orientation and motion. 
Parallel motion gave a single continuous solution curve. At some transition point between 
the parallel and orthogonal conditions, the solution breaks into two distinct non- 
continuous branches. The transition point is a singularity and makes it difficult to 
characterize algorithm properties such as the number of required sensors to guarantee 
unambiguous geolocation solutions. Furthermore, the symmetric properties of the FDOA 
solution vary. In the parallel motion condition, we found a symmetric solution about the 


sensor axis. With orthogonal motion, this symmetry fails. 
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Figure 16. | Two FDOA sensors geolocate in the orthogonal motion condition. 
3. 60° Sensor Motion Condition 


We consider one case in the aforementioned transitional region between parallel 
and orthogonal sensor motion relative to the sensor axis. The sensors and emitter are 


located at the same positions described previously. Let us choose, however, the sensor 
velocity v, =v, =(5.00%+8.663) m/s. Such a velocity forms a 60° interior angle with 
the sensor axis (+X orientation). The resulting constant-FDOA solution curve is given in 
Figure 17. The solution is again a non-continuous, multiple-branch curve. The angular 
offset effectively skews the direction of the two lobes of the curve. The extent of the 


FDOA solution dependency on sensor motion direction is readily apparent. 
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Figure 17. | FDOA sensors geolocate for 60° motion relative to the sensor axis. 


E. THREE-DIMENSIONAL SOLUTION 


The mathematical solution for the FDOA equation in three dimensions follows 
closely that of the two-dimensional case. While not increasing in theoretical complexity, 
the parametric nature of the proposed solution raises the computational demands 
significantly. In two dimensions, the x-y plane is described by the fixed polar angle 


@=90°; consequently, the previous solution only demanded computation time on the 


order of O (m) , where m is the number of azimuthal angle samples in the range [0, 2r| 
and O( -) represents the traditional big O notation. Extension into three-dimensional 
spherical coordinates requires n polar angle samples in the range [0°,180°] per azimuth 


angle; therefore, the associated computation proceeds as O(mn). For any reasonable 


degree of spatial precision, significant computational resources are required. 


Accordingly, it is increasingly difficult to render high-quality three-dimensional images 
af 


of the FDOA solution surface. Thus, in the figures that follow, only discrete sample 
points are plotted. An effort has been made to choose a sufficiently small spatial sample 


size so as to capture the general shape of the FDOA solution surfaces. 


1. Parallel Motion Condition 


Two sensors 5, and 5, are located at coordinates (—5,0,0) and (5,0,0) km and 
have identical velocities of +10 m/s. A stationary emitter is placed at coordinates 
(6, 4,0) km. The corresponding FDOA solution surface is presented in Figure 18. On the 


left, a three-dimensional view is given. On the right, a top-down view is given from the 
perspective of the +z axis; all points are therefore projected onto the x-y plane. The 


sensors and emitter are depicted as red and black markers, respectively. 
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Figure 18. | Three- and two-dimensional views are given of the sampled FDOA surface 
for the parallel motion condition. 


Here, the symmetry about the sensor axis is convenient when generating the 
three-dimensional surface. Indeed, the solution is simply a surface of revolution using the 
two-dimensional “peanut’-shaped solution found previously. The axis of revolution is the 


sensor axis; here again, since both sensors are located directly on the x-axis, the sensor 
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axis is the x-axis. This is the simplest condition for which the FDOA equation may be 
solved in three dimensions. 


2s Orthogonal Motion Condition 

Let us reset the sensor velocities in the scenario above to +109 m/s. The sensors 
now move orthogonally with respect to the sensor axis (x-axis). With all other 
environmental parameters held constant, the corresponding FDOA solution surface is 
given in Figure 19. On the left, a three-dimensional view is given. On the right, a top- 
down projection onto the x-y plane is given. Now the axial symmetry is lost. The solution 


surface is no longer a surface of revolution, and the need for a parametric solution 


strategy is validated. 








16 
16 
5 
a5 2g & 
5 6 & y 
a OO “a 0 e 
: @ : 
us A 
-5 
-10 16 
16 0 
0 | -10 
y-axis (km) -10 -10 x-axis(km) -10 0 10 
x-axis (km) 


Figure 19. Three- and two-dimensional views are given of the sampled FDOA surface 
for the orthogonal motion condition. 
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3. 60° Sensor Motion Condition 


Lastly, we revisit the 60° sensor motion condition, where the sensor velocities 
are set to (5.00% + 8.665) m/s. The FDOA surface is shown in Figure 20. As before, a 


three-dimensional solution is given on the left, and a two-dimensional projection in the x- 
y plane is given on the right. Clearly, no symmetry about the sensor axis (x-axis) is 
available; thus, the surface is not the result of revolution. Although the two lobes found in 
the solution are similar in shape, the symmetry is more difficult to define mathematically 


and offers little advantage in easing the solution requirements. 
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Figure 20.‘ Three- and two-dimensional views are given of the sampled FDOA surface 


for the 60° relative motion condition. 
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V. NOISE CONSIDERATIONS 


A. SOURCES 


As with any real system, noise is an unavoidable complication. Here, the term 
noise refers to any mechanism by which errors are introduced into a system. Some 
geolocation systems may be particularly tolerant to various types of noise, yet the 
corruptive impact of noise is felt to some degree in all facets of a geolocation algorithm. 
One might broadly categorize noise sources as internal or external to the geolocation 
system. External errors may arise from such causes as inhomogeneity in the propagation 
medium, thermal agitation of the medium, nonzero dispersion characteristics of the 
medium, and nonlinear delay profiles for LOS transmissions. For satellite systems, with 
significant propagation distances, issues such as scattering and refraction in the 
ionosphere greatly impact signal quality and delay. Internal errors may arise from thermal 
recelver noise, inaccurate synchronization to a global clock, assumption inaccuracies, 
hardware degradation and failure, or computational rounding error. Altogether, the list of 
potential sources of error is extensive. Thus, accurate statistical characterization of such a 


complicated system through a simple noise model is a challenge, to say the least. 


B. NOISE MODELS 


With so many diverse noise sources and environmental parameters to be 
considered, developing a representative geolocation noise model is an area of significant 
attention. For instance, the Community Coordinated Modeling Center (CCMC)—a multi- 
agency space weather partnership led by NASA—provides public access to nine distinct 
ionospheric computational electromagnetic (CEM) models developed by a variety of 
universities and research institutions [18]. Modeling these propagation characteristics 
alone can be extremely complicated, so in order to maintain our research focus of data 


fusion in time, we initially choose a basic additive noise model. 
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1. TOA Noise Model 


For time-domain measurements, we choose to model the net environmental 


impact and system errors using a single exponentially-distributed, additive noise 
parameter specific to a given sensor s, . The exponential distribution is selected to model 
the time delay for two reasons. First, the exponential distribution places a hard lower 
bound on the delay; the beacon signal cannot arrive at sensor s; with delay smaller than 
the propagation time associated with the free-space propagation speed c,. Second, the 
exponential distribution is relatively simple to handle mathematically, which helps avoid 


unnecessary complications in the analytical work to follow. 


To clarify the notation convention used hereafter, an exponential random variable 


x with the probability density function (PDF) 


7 0, x<0 55 
ie (x)= Ae x>0 ( ) 


given in terms of the statistical parameter 2 may be alternatively defined using the 


condensed notation x ~ Exp(A). This generic random variable x has mean 4, =A", 


variance 0,7 =”, and standard deviation 0, =A" [19]. 


Let us define the time-domain noise parameter associated with sensor s, as 
6,, ~ Exp(A,). A TOA measurement collected by sensor s; may then be expressed as the 
sum of the theoretical deterministic time delay (see Equation 1) and the random variable 
0,;: 


deterministic 


random 
mi) 
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The resulting noisy TOA measurement is a random variable with mean 











By Aiea (57) 
H, = Cy wy 


and standard deviation 
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C= (58) 


Indeed, the exponential distribution is a rough first approximation of the time- 
domain noise effects on a TOA beacon signal. Still, it provides a suitable first 
approximation of the stochastic properties associated with atmospheric propagation and 
noise without introducing excessive complexity into the problem. In light of this 
assumption, we emphasize the additive random nature of the noise parameter used in this 
model. This modular design easily allows for future modifications. If another PDF— 
defined either analytically or empirically—is found to better suit the propagation model, 
one need simply replace the PDF description and re-execute the simulations included in 
this research. Alternative distributions that have found applications in the wireless 
environment are the Rayleigh, Rician, and lognormal distributions [20]. Future work 


might explore the incorporation of these distributions in this geolocation model. 


2. TDOA Noise Model 

The TDOA noise model is largely built upon the TOA noise model. The TDOA 
measurement is defined as the difference between the two distinct arrival times at sensors 
s, and s, (see Equation 3). The arrival times 7, and t,may be decomposed into one 
deterministic and one random-valued term, as in (56) with the TOA model. We thus 
express the TDOA measurement as 


deterministic 4 
random 
& J | | 


r= (=" “3 ae 5 =—(b |-f.-7[))+(8:-6,). 69) 
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Let us assume the random noise parameters 6,, and 6,, are independent and identically 


distributed (IID) exponential random variables, implying 4, =/2,. When computing the 


difference between two independent exponential random variables, we subtract the 
respective means and add the respective variances to characterize the resulting 


distribution. If the random variables are further identically distributed, the resulting 
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distribution has mean yz, —/4,=0 and variance o,’ +0," =20,”. Applying this to (59), 


we find that the TDOA measurement Tj is arandom variable with mean 











I-[r.-7l) @) 


1 
bei =. —( r, TT, 
Co 


and standard deviation 


C,,, =. (61) 


Note how the same time delay distribution characterized by the parameter A, causes a 


factor of 2 increase in the TDOA measurement standard deviation compared against 


the TOA measurement standard deviation. 


3. FDOA Noise Model 


Since the FDOA algorithm is built upon variations in the Doppler shift, we must 
consider the effect of environmental and system noise on the Doppler shift of a given 
signal. In the frequency domain, a propagating signal does not have the same restrictions 
as in the time delay model, where there is a hard minimum delay associated with the free- 
space propagation speed. Rather, it is equally likely that the observed frequency of arrival 
at a given sensor will be less than or greater than the anticipated center frequency. As 
such, we choose to model the Doppler shift noise as a Gaussian-distributed random 
variable, which is symmetric about its mean value. To be explicit, a generic Gaussian 


random variable x, with mean yz and standard deviation O, is characterized with the 


PDF [19] 


(x-u) 
ae (62) 
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We introduce the condensed notation which alternatively defines the same Gaussian 


random variable x ~ N ( HM, o*) : 


Recall that the FDOA measurement is the difference in Doppler shift experienced 


by the emitter beacon signal upon arrival at sensor s; and s, (see Equation 36). Let us 
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model the Doppler shift observed at sensor s, as the sum of the theoretical Doppler shift 
due to the radial component of the relative velocity of the emitter (see Equation 32) and a 


. . . . 2 - 
Gaussian-distributed random variable 6, ~ N ( sao F op ) ; 


Af, = fyi )+5p- (63) 


Co 


We similarly model the Doppler shift observed at sensor s, in terms of the Gaussian- 


distributed random variable 6, ~ N (Us,+0% rs 


Ne fe aa vi, )+6 


fi* 
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(64) 


The FDOA measurement Af;,=Af;-Af, computed between sensors s, and s, may 
now be expressed as the sum of the deterministic FDOA value predicted by (41) and a 
random noise component: 


deterministic 
random 


Jo (yi = 65 
ae (vi Ws i \+ (6; 5). (65) 
0 
The random noise component of the FDOA measurement, given by (6 Ba Oe Is is 
the difference of two Gaussian distributions; to characterize the resulting distribution, we 
subtract the means and add the variances. The resulting random variable is Gaussian- 
distributed as (6 RO s a N( Msq— Us goFs es +0; be ) . Considering the deterministic 
FDOA value as a bias, the total FDOA measurement is Gaussian-distributed according to 
fi : VE vi + sq — Ms 5 Oa: +055 *) . (66) 
0 
This is a generalized form of the FDOA noise model. We assume the Doppler shifts Af, 


and Af, are IID and zero-mean. In this case, Cres =0; ef =o,, , and (66) reduces to the 


uP 


more manageable form 
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Totes 2 
Af, -[ 2-4). 24 (67) 
Notice that the standard deviation of the FDOA measurement is scaled by the factor /2 
when compared with the Doppler shift observed at an individual sensor. Like in TDOA, 
an FDOA measurement is the difference of two independent random variables; as such 


the total variance is larger than that of either individual random variable. 


C. NOISE PARAMETER ESTIMATION 


To complete the TOA, TDOA, and FDOA noise models, parameter values are 
needed for the time delay and Doppler shift noise profiles. Two broad cases are 
considered. First, the sensors are modeled as high-altitude geolocation satellites. Second, 
the sensors are modeled as low-altitude UAVs. Depending on the altitude, the time delay 
and Doppler shift profiles are significantly different. 


1. High-altitude Sensors 


A significant source of signal deterioration in high-altitude satellite 
communications is propagation through the ionosphere. The ionosphere is a series of 
layers in the Earth’s upper atmosphere with characteristically high concentrations of 
ionized particles and free electrons. The altitude boundaries of the ionosphere are highly 
variable, depending on solar weather conditions and the time of day. Typically, the 
ionospheric layer ranges from a lower altitude of 50-600 km to an upper altitude of 600- 
2,000 km [21]. All altitudes are referenced to the Earth’s surface unless otherwise stated. 
An electromagnetic (EM) signal passing through this thick ionized medium may 
experience significant degradation from propagation effects including attenuation, 
Faraday rotation, group delay, scintillations, and refraction [22]. The task here is to 
determine the aggregate impact of such ionospheric effects on the time delay and Doppler 


shift of a transiting beacon signal. 


The following discussion uses the orbital and communications figures attending to 
the Global Navigation Satellite System (GNSS) satellite constellation. The orbital altitude 


is 19,100 km, the orbital speed is 3.95x10° m/s, and the primary communications 
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operating frequency is within the L1 band at 1575.42 MHz [23], [24]. The emitter is 


assumed to be stationary. 


a. Time Delay Variation Characterization 


Let us begin our time delay variation discussion with a brief review of EM wave 
propagation theory. The following is adapted from the extensive treatment given to 
Maxwell’s equations and plane-wave propagation in [25]. As with any electromagnetic 
field, Maxwell’s equations govern a satellite communication signal transiting the 
ionosphere. The homogeneous wave equation is a solution to the time-harmonic 
Maxwell’s equations and gives 

VE+o uk =0, (68) 
where @ is the angular frequency of the oscillating field; 4 and ¢ are the magnetic 
permeability and electric permittivity, respectively, of the medium; and E is the 
monotonic electric field vector in phasor form. The wave equation is frequently 


expressed in terms of the wavenumber 

k =a ue. (69) 
Often, the permeability is written in terms of the relative permeability “4, and the free- 
space permeability 4, : 


H= Ly. (70) 


The permittivity is likewise written in terms of a relative and free-space value as 
E=E,Ey. (71) 


We may also define the index of refraction n such that 


n = H€,. (72) 
Except for very rare types of materials, n > 0 and the positive root of (72) is used. Noting 
that the speed of light is defined as c, = / a LoEq » the phase velocity of the electric field 


is 
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y= (73) 


This is the velocity of the constant-phase points of the traveling EM wave and may 
exceed the speed of light in a vacuum c, in some instances without violating physical 
laws. 

When an information signal is modulated onto an EM carrier, the result is no 
longer monotonic, and the bandwidth is nonzero. Consequently, the information 


propagates through the medium at what is termed the group velocity v,, which is not 


necessarily equivalent to the phase velocity v, and is bounded by c, as a maximum 


value. The net result is a delay in the arrival of the information signal that is termed the 
group delay Az,. The analytical solution for the group delay involves the use of the 
Appleton-Hartree equation and is too involved for the scope of our present discussion; 
see the chapters entitled “Plane Waves” and “Ionospheric Propagation” in [22] for an 
extensive treatment of the group delay. Instead, we use the final approximate result for 


group delay when propagating through an ionized medium: 
1 oe 
Ag ee ie (74) 
82° \ Em,c, ) f 


where e is the fundamental charge of an electron, m, is the mass of an electron, f is the 





frequency of the EM wave, and T.. is the total electron content of the transmission path 


[22]. We use the average group delay as the statistical mean of the time delay noise 


profile in our noise model, but first we must solve for all variables in (74). 


The approximate group delay result above is based upon a model of the 
ionosphere as a Lorentz plasma with volumetric electron density N, [22], [26]. The total 


electron content is defined as the path integral of the electron density: 


T,, =| N, ds. (75) 


Let us assume the electron density is given by the simple piecewise expression 
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0, h<h. 
N,= 5N h <h<h- 


e ion? ion ion ? 


0, h>h* 


ton 


(76) 


where N, is the constant electron density within the ionosphere and h, and h, are 


ton ton ton 


the lower and upper altitude limits of the ionosphere. If we further assume vertical 


propagation from an emitter entirely below the ionosphere (h < hin) to a satellite beyond 


the ionosphere (h Se ) , the total electron density reduces to 


ton 


Toe = Nin (Tibn — Hon) (77) 


Recall from above that the lower altitude limit of the ionosphere typically ranges 


from 50-600 km; let us choose 325 km as the average value of h,, 


ton 


. The upper altitude 


limit ranges from 600—2,000 km; choose 1,300 km as the average value of h;. We thus 


ton 


say, on average, the thickness of the ionosphere is (h; — Hin =975 km. 


ton 


We determine an electron density of the ionosphere from computer prediction 
models. One model of interest is the Global Assimilation of Ionospheric Measurements 
(GAIM) developed at Utah State University. Recent simulation run results are made 
available to the public by the CCMC [18]. Consider a typical simulation forecast run for 
January 2, 2012. For an altitude of 844 km, the electron density N, is forecasted at 15- 


minute intervals for all latitude and longitude coordinates. The simulation results for 
00:00:00 UTC are shown in Figure 21. Here, the electron density is plotted in 


electrons/em® with a colorbar on the right side of the figure for scaling purposes. 
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Figure 21. The electron density of the ionosphere at altitude 844 km varies with latitude 
and longitude on January 2, 2012 at 00:00:00 UTC. From [27]. 


Results from 12 hours later at 12:00:00 UTC are shown in Figure 22. Again, the 


electron density is plotted in electrons/cm’ and a colorbar is provided. 
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Figure 22. The electron density of the ionosphere at altitude 844 km varies with latitude 
and longitude on January 2, 2012 at 12:00:00 UTC. From [27]. 
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Notice the severe dependency of N, on the time of day. As the Earth rotates a 


half-revolution in the 12-hour interval between these two snapshots, those longitudes 
with areas of high electron density (red) shift to areas of low electron density (blue), and 


vice versa. Note further that high electron density regions are mostly grouped around 
equatorial latitudes. Since there are such large variations in N,, we design for the worst- 


case scenario. Here, the highest electron density induces the longest group delay, so we 


select the highest measured electron density in either sample to characterize the 


ionosphere, namely N, =172500 cm™® =1.725x10'' m®. 


Under these basic assumptions, the average total electron content of the 


ionosphere is evaluated from (77) as 


T,, = (1.725x10"' m™*)(975x10° m)=1.682x10"7 m™. (78) 


We now evaluate (74) for the average group delay due to ionospheric propagation at the 


GPS L1-band operating frequency of 1575.42 MHz: 


1 Oey 
At, ~1(—|& 
1 (1.602x10-? cy 
81° | (8.854x10°" F/m)(9.109x10! kg)(2.998x10* m/s) ‘ (79) 








(1.682x10" m”~) 
(1.5754210" Hz) 





=9.110 ns. 


This average group delay is used to characterize the exponentially-distributed time delay 


noise 6,, ~ Exp(A =(9.110 ns) '} 


We now complete the time-domain noise model for high-altitude satellites. The 


noisy TOA equation given in (56) may be expressed as 


deterministic 
random 





T, = Ets Exp(a=(9.110 ns) '). (80) 


i 
0 
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Similarly, the noisy TDOA equation given in (59) may be expressed as 


deterministic random 


|r. —r,l))+(Exp(2 =(9.110 ns) ')- Exp(2 =(9.110 ns) ')), (81) 














rT 








1 
Tt, =— 
b. Doppler Variation Characterization 


The ionosphere also impacts the Doppler shift of a beacon signal via two primary 
mechanisms: refraction and phase alterations due to the complex refractive index of the 
ionosphere as detailed by the Appleton-Hartree equation. In fact, the frequency shift 


exerted by the ionosphere on a signal may be expressed in similar terms as the 


previously-discussed group delay Az, (see Equation 74) [28]: 


1 e 1 \d 
ras = Sr (<4) (7). (82) 


The time-derivative of the total electron content implies that the ionospheric Doppler 





shift is nonzero only when 7). is a time-varying quantity. We also note that Af, is 


inversely proportional to the operating frequency /f,. Thus, ideally, higher operating 
frequencies will suffer less from ionospheric Doppler shifts. 

Willman has collected experimental data on the ionospheric Doppler shift due 
to refraction. For a satellite signal transiting the ionosphere at an operating frequency 
of 54 MHz, he recorded a maximum Doppler shift of 40 Hz [29]. This is the worst- 


case scenario he recorded. We scale this error to the GPS operating frequency of 


1575.42 MHz according to the inverse-frequency dependency of (82): 


af, =| 20 |a ft, =| 22 ME2 _ |(40 wz) =1.371 Bz. (83) 
Te 1575.42 MHz 
Let us assume this ionospheric Doppler shift is the standard deviation of the zero-mean 


Doppler shift noise variables 6,, in (63) and 6, in (64). 


We complete the frequency-domain noise model for high-altitude satellite 


sensors. The noisy FDOA equation given in (65) may be expressed as 
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deterministic random 


af, =2(vi,-v))+ (1v(0.@.371 Hz)’ )-(0,(1.371 Hz)’)). (84) 


Co 





The noisy FDOA measurement may be alternatively expressed in condensed notation as 


Af, ~ w(t -vi,), 2(1.371 Hz | (85) 
0 
2. Low-altitude Sensors 


An alternative to deploying a permanent satellite constellation is using mobile 
UAVs or other airborne platforms as geolocation sensors. While such systems do not 
have the longevity or coverage area of a satellite constellation, low-altitude sensors have 
the advantage in flexibility and maneuverability. If the ionosphere was the main 
propagation concern for GHz-range satellite links, multipath fading is the primary 
mechanism of signal degradation for low-altitude sensors. As with all wireless 
communications links, the unique characteristics of the environment have a severe impact 
on signal propagation. Unlike ionospheric distortion, which was highly variable yet 
relatively predictable using forecast models like GAIM, multipath fading is highly 
variable and very difficult to predict. Significant research has been directed toward 
modeling multipath effects on delay spread and Doppler spread. Analytical solutions are 
far too complicated for practical usage. Instead, numerical computations or empirical 
approximations are typically used to estimate multipath effects [20]. Here, we turn to the 
literature to gather reasonable multipath parameters for low-altitude UAV 


communications links. 
For the following discussion, we assume all sensors are UAVs at an altitude of 
10,000 ft (3,048 m) and traveling at a cruising speed of 100 knots (51.44 m/s). All 


pertinent sensor communications operate within the unlicensed industrial, scientific, and 


medical (ISM) band at 2.4 GHz. The emitter is assumed stationary. 


a. Time Delay Variation Characterization 


As a signal propagates radially from a transmitter, some portion of the signal 
power propagates along the LOS path and some is lost due to spherical spreading of the 
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wave in the medium. When the environment is filled with objects to reflect, scatter, and 
diffract this wayward signal power, some may be redirected back toward the intended 
receiver. Upon arrival, there will be an associated time delay due to the circuitous path 
taken en route. The resulting difference in propagation delay between the LOS signal and 
the multipath signal is termed the delay spread. Let us define the root-mean-square 


(RMS) delay spread as [30] 





DA ine). 
YB, 


where 7, is the delay spread, f is the signal power versus delay profile, and 7, is the 


AT 


Ss 





(86) 


mean delay spread value defined as [30] 


7, = (87) 


Newhall et al. collected experimental data for delay spread of a 2 GHz signal 
versus elevation angle above the horizon. They found that delay spread falls off 
significantly as the elevation angle increases, suggesting a large dependence on the high 
concentration of particulate matter in the atmosphere at low altitudes. At a maximum 
elevation angle of 30° above the horizon, they observed an average RMS delay spread of 
18.3 ns [30]. For this model, we shall exclude such low elevation angles for the sake of 
accuracy. At higher elevation angles—nearer to vertical-path propagation than 30°—the 
delay spread is assumed to further decrease. Accordingly, we choose an average RMS 


delay spread of At, =10 ns to characterize the time domain signal delay variation. 
For low altitude UAVs, the noisy TOA equation given in (56) may then be 
expressed as 


deterministic 
random 

















r, = Hl + Exp(4 =(10 ns) '), (88) 


i 
0 
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The noisy TDOA equation given in (59) may be similarly expressed as 


deterministic random 


r,—r,| —lr. -r,|) +(Exp(2 =(10 ns)" - Exp(4 =(10 ns)" )). (89) 














1 
T=— 
=H 
b. Doppler Variation Characterization 
Consider a UAV with onboard accelerometers and inertial navigation equipment 
that has a speed estimation tolerance of +0.1 % . At a cruising speed of 100 knots (51.44 


m/s), this corresponds to a speed uncertainty of +0.0514 m/s. If we assume the worst- 


case scenario where the sensor is moving directly toward or away from the stationary 


emitter, (32) gives an average Doppler variation due to speed uncertainty of 





9 
Af, = L(y) = (ae) +0.0514 m/s) =+0.412 Hz. (90) 


Cy (2.998 10" mis) 
As discussed previously, we assume the Doppler variation about the deterministic mean 
value to be Gaussian-distributed. Let us use this Doppler variation as the standard 


deviation of the zero-mean ITD random variables 6 and O f° 


We may thus rewrite the noisy FDOA equation given in (65) as 


deterministic random 





Af, = Av -vi,) +(1v(0,(0.412 Hz)")—(0,(0.412 Hz)')), (91) 


If the deterministic FDOA component is considered a constant bias, the total FDOA 


measurement may also be expressed in condensed notation as 


Af, ~ n( 200 -vi), 2(0.412 H)| (92) 


Co 


a) 
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VI. SYNTHETIC APERTURE: TOA 


A. OVERVIEW 


The principle of a synthetic aperture is straightforward. A series of measurements 
or data samples in time are compiled together in such a way as to yield a refined high- 
resolution measurement. The implementation of the synthetic aperture concept is often 
more nuanced. Significant research in the past fifty years has been devoted to using 
synthetic apertures in various fields. Airborne radar systems implement SAR techniques 
as they fly over a target of interest [31]. Synthetic aperture sonar (SAS) is used in 
maritime environments to generate high-quality underwater imagery. As a vessel travels 
along the ocean, it collects sonar returns at either regular or irregular intervals. These 
returns are subsequently synthesized using a SAS algorithm to provide a refined image of 
the underwater environment [32]. The synthetic aperture concept has made inroads as 
well with the medical ultrasound imaging field. Synthetic aperture ultrasound images are 
refined enough to capture blood flow within an organism [33]. We aim to replicate the 


synthetic aperture principle in geolocation algorithms. 


B. SOLUTION STRATEGY 


Recall from Chapter II that two TOA sensors in three-dimensional space can 
localize an emitter position to a circular curve, namely the intersection curve between the 
two TOA spheres (see Figure 5). The TOA algorithm is unique in this sense; for any 
spatial sensor orientation, the circular solution curve may always be defined analytically. 
We exploit this fact in the following heuristic solution strategy for the synthetic aperture 


TOA problem. 


1. Parametric Circular Solutions 


An analytically-defined curve may be efficiently adapted to a discrete numerical 
computation environment through parametric expression. As discussed above, at any 


given time step, two TOA sensors can localize an emitter to a circular curve in three- 


a7: 


dimensions. Thus, one must know the parametric equation governing the circular emitter 
curve. 
A circle with center coordinates given by the position vector r, and radius p 
may be expressed in terms of the parameter ft by the three-dimensional position vector 
p(t)=r,+ pcos(t)u+ psin(t)y, (93) 
where # and ¥ are arbitrary orthogonal unit vectors parallel with the plane in which the 
circle lies [34]. Let this plane be characterized by the unit normal vector k . It then 
follows that 
k=+(ax?), (94) 
where the (-)x(-) operator denotes the cross product. Note that (93) is periodic with 


respect to parameter ¢. Thus, we can fully define the circle while restricting f¢ to the 


domain [0°,360°). A generic diagram of such a parameterized circle is given in Figure 


ods 





Figure 23. A generic circle in three dimensions is shown as the parametrically-defined 
position vector p(rt). 
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The next task is to generally express the two-sensor TOA solution in its 


parametric form. Accordingly, the circle radius, center coordinates, and unit vectors are 


required. Consider the following noise-free scenario. Two TOA sensors s, and s, are 
located at the three-dimensional coordinates specified by the position vectors r, and r,, 


respectively. The unit vector F,, points in the direction from sensor s, to sensor s,. An 
emitter e generates a beacon signal while located at the coordinates specified by the 


position vector r,. The sensors receive the TOA signal and calculate the corresponding 


TOA spheres. From the TOA measurements, the scalar range values Ir. and 














r,,|| are 


defined according to Equation 1. A cross-section of this three-dimensional scenario is 
taken such that the cross-sectional plane contains both sensors and the emitter. Such a 


plane is shown in Figure 24. 


sensor 
axis 





Figure 24. A cross-sectional plane of a two-sensor TOA solution is shown. 
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The circular solution is the path traced by the point e when revolved about the 


sensor axis. Therefore, denotes the radius of the circular solution, and the point c 


denotes the center of the solution circle. To express these parameters in terms of known 


quantities, we first apply the law of cosines: 


2 2 2 
lee =r, -7 + = 




















al 














al 








lr, -1,||cos(,). (95) 
Bounded on the range [0°,180°], the angle @, may be expressed as 


2 2 





















































2 
r,—r,|) +r, -(r 
0, = arccos ro=nll +H 2 (96) 
2 frills 
The circular radius is then given as 

p=lr.,||sin(@), (97) 

and the position vector of the circle center is given as 
re =r, +| Py cos(@,) |F,,. (98) 














To define the required parameterization unit vectors, we first observe that the 


solution circle lies in the plane perpendicular to the sensor axis. Therefore, we choose 


k=f,. (99) 
We now must choose two additional unit vectors W@ and ¥ that lie within the plane of the 
circle and are mutually orthogonal with k . There is an infinitude of such combinations, 


so we arbitrarily choose #@ such that it has a z-component of zero. This implies that 


AK 


fi-P,, =(1,£+0,9+08)-(A,F+A,,9+%,2)=0. (100) 
It may be shown (see Appendix) that one of two real roots for & is given as 


a 2 na 2 


L r 
ipl 21 A A A 
u= aD) > 5 x+ {tS j+(0)Z. (101) 
iy +h), Nix +h, 


Lastly, the unit vector ¥ may be defined as 


p=uxk, (102) 


where k and w are given in (99) and (101) in terms of known quantities. We can now 
fully express the parameterization of a circular TOA solution curve in the form of (93) 


using known, assumed, or measured quantities. 


Recall that the discussion so far has pertained only to one single instance in time. 


If the position vector p, describes the circular solution of two TOA measurements at a 
given time step n and no noise affects the system, then there is some parameter value ¢ 
for which p,, (z' ) accurately specifies the emitter position. To determine the specific 


value of f,, we shall use the subsequent circular solution at time step n+1. 


2. Circular Intersections 


The synthetic aperture concept presupposes that data from two or more distinct 
instances in time may be compiled in such a way as to refine the overall solution. In the 
context of the TOA algorithm, this involves finding the intersection of two or more 
circles in three-dimensional space. With the known parametric equations of subsequent 


circular emitter curves in time, this is simply a matter of solving for the theoretical 


r 


parameter values of f, and f,,, for which p, (t') = p,,, (¢/,,). 


The process is more complicated when noise is introduced into the system. Now, 
subsequent noisy circular emitter curves do not necessarily contain the exact coordinates 
of the emitter and do not necessarily intersect one another. We circumvent this issue in 


our algorithm by finding the nearest two points on the two circles. This is accomplished 


using numerical computation; we perform an exhaustive search through all values of ¢, 


and ¢ 


n+l 


to find the parameter pair (t”,7",,) such that 


| Py (ty) Pyas (toss) |: (103) 


Note that we designate this parameter pair with double-primed notation to emphasize that 


" 


P, (17) 7 Pasi (a )| = min| | 








the associated points do not necessarily intersect. There is some inherent quantization 


error introduced into the problem since this search is implemented numerically; however, 


if a sufficiently large number of samples is used to discretize the parameters ¢, and ¢ 


n+l? 
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then the quantization error becomes negligible compared with the effect of the noise due 


to time delay of the TOA beacon signal. 


To compute the estimated position of the emitter, we now find the midpoint of the 
two nearest points on the two circles. Let us represent the estimated emitter coordinates 


based on measurements at time steps n and n+1 with the position vector r”,, defined as 
p p é 


yf? 


roy = Pell Poult) (104) 


For a synthetic aperture composed of N time steps, the refined emitter location estimate 


r; is computed as the spatial centroid of all N-—1 available r/, estimates. 


Mathematically, we express the refined estimate as 
1 N-l 
r,=—— ) r. (105) 
N a n=1 i 


This refined emitter estimate is the result of the synthetic aperture TOA algorithm. 


3. RMS Error 


To assess the performance of the algorithm under various sensor orientations and 
synthetic aperture sizes, we must quantify the associated error performance for a given 
noise profile. The distance error of a single synthetic aperture estimate relative to the true 


emitter position is given as 


re 














" 
rT? 


(106) 
To characterize the general noise performance, it is customary to consider the 


theoretical MSE of our estimator algorithm, which is defined as 














MSE = E[d? |=2| |r"-r, a (107) 
where E|-]| denotes the expected value operator [35]. The MSE here carries the units of 


square meters and has limited physical significance. A quantity of more relevance is the 


root-mean-square error (RMSE) and is simply given by 


RMSE = VMSE = || d’ |. (108) 
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The RMSE carries units of meters and is more useful for physical reference. 


Ideally, the RMSE of an estimator algorithm is determined analytically, and the 
expected value operation is evaluated in closed form. For complicated estimators, this is 
increasingly difficult. Instead, the expected value operation is frequently approximated by 
the average squared error over many independent trials. In our numerical computation, 


we assume the approximation 


RMSE # ,/—)d,,” (109) 
M 
holds for sufficiently large M . 


Cc. RESULTS 


Consider the following scenario. A stationary emitter e is located at the 


coordinate system origin and generates a TOA beacon signal. Two TOA sensors 5s, and 
S, are initially located at the Cartesian coordinates (—b- Ab,-b,h) and (-b,-b,h), 
respectively. Each sensor travels at the same speed v,; the direction of the sensor 


velocities v, and v, may vary. Over time, the two sensors collect N total noisy TOA 
measurements for inclusion in the synthetic aperture estimate—we refer to this parameter 
as the size of the synthetic aperture. The time delay standard deviation o-,,, is 
incorporated into the TOA measurement as detailed in Chapter V. The constant time 


interval between measurements AZ;, is chosen such that, at the end of the total synthetic 


aperture time duration T,,=(N-1)AT,,, the sensors have traveled distance 2b. This 


SA ? 


necessitates that 


_v,(N-1)AT, 


b= : (110) 


Establishing a generalized setup geometry such as this allows us to maintain some 


comparability between the following simulations run against different system variables. 
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Three primary variables are considered in this discussion: the angular spread 
between sensor velocities v, and v,, the size N of the synthetic aperture, and the 


altitude h of the geolocation sensors. Two distinct cases for sensor altitude—satellites at 
altitude 19,100 km and of UAVs at altitude 10,000 ft (3,048 m)—are simultaneously 
presented in the following two sections regarding sensor velocity angular spread and 


synthetic aperture size. 


if Sensor Velocity Angular Spread 


The sensor velocity angular spread refers to the angle @ between the sensor 


velocities ¥, and v,. It may be defined in terms of a dot product as 


oatcos| Hi Jasco PP) (111) 
Iv: lv v, 


if a is confined to the range [0°,180°] . Thus, to sweep through a range of spread angles, 





we must vary one or both of the sensor velocity directions. 


Let us choose a fixed sensor velocity v,=v,y. To vary the spread angle, we 

define the variable sensor velocity 
v, =v, sin(a)X+v, cos(a@) ¥, (112) 
with @ set to vary on the range [0°,90°]. A top-down view of this scenario is diagramed 


for illustrative purposes in Figure 25. 
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Figure 25. — The generic sensor velocity angular spread scenario is diagramed. 


The TOA circular emitter curve generated from each pair of TOA measurements 
collected at a given time step is discretized using 20,000 sample points. To further boost 
the resolution of the estimated position and improve on computation time, we impose a 


maximum altitude h on the emitter search domain. This constraint effectively 


condenses the 20,000 sample points on the TOA circle to the region of interest rather than 


wasting computation time by considering points known to be very far from the true 
solution. If the emitter is known to be at high altitudes, h,,,, and the number of sample 


points may be increased at the expense of computational efficiency. The principle is 
illustrated in Figure 26; black dots represent discretized sample points on the TOA circle, 


and the red dot is the true emitter position. 
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TOA circle 





Figure 26. An altitude constraint is imposed on the TOA emitter search, condensing the 
search to an area of interest. 


The synthetic aperture time intervals AT,, were chosen according to the following 
process. A maximum total synthetic aperture duration time Jj, for satellite sensors was 


set to 15 minutes. With N =5 aperture samples, we find that AT), =225 seconds for 


satellite sensors. To ensure that the system geometry is comparable for UAV sensors, the 
time interval was scaled according to sensor altitude and velocity to keep the swath angle 


of the sensors constant (see Appendix for details on scaling and swath angle). The result 
is that, for N=5 aperture samples, AZ;, is 2.72 seconds for UAV sensors. Other 


pertinent simulation parameters have been collected and organized below in Table 1. 
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Table 1. Angular spread simulation parameters are given for the 
synthetic aperture TOA algorithm. 


Satellite Sensors UAV Sensors 
19.10 10° 3.048 x 10° 








3.889 <10" 51.44 





a (degrees) variable on [0,90] variable on [0,90] 





N (samples) 5 5 





225 2.72 





1.750x10° 280.1 





1000 100 




















Each simulation was run under two noise conditions. A high-noise environment 
was simulated using the time delay standard deviation o,,, listed in Table 1. To model a 


low-noise environment, the standard deviation is reduced to 10% of the listed value. This 
can model a nighttime environment with minimal solar agitation and disturbance in the 
ionosphere for the satellite scenario. For the UAVs, at nighttime, thermal convection 
currents in the lower atmosphere may be less pronounced and, thus, have a smaller 
impact on the delay of a radio frequency (RF) signal. In all cases, we iterate the algorithm 
ten times to help generate an approximation of the RMSE. More iterations are preferable, 


but the significant computation time of the algorithm made this prohibitive. 


a. Satellite Sensors 


The RMSE for the synthetic aperture TOA algorithm applied to satellite sensors is 
plotted against sensor velocity spread angle @ in Figure 27. All data points are plotted as 
circles centered on the data value, and all traces between data points are generated using 
basic linear interpolation. We do not apply any regression analysis since the sample size 


is so limited. 
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Figure 27. | RMSE is plotted against sensor velocity angular spread q@ for satellite 
sensors. Dashed traces are linear interpolations between adjacent data points. 


These satellite sensor simulation results highlight several important characteristics 
of the synthetic aperture TOA algorithm. Note the order of magnitude variation in the 
high-noise trials as the sensor velocity spread angle is varied. Also, observe the general 
trend of decreasing RMSE with increasing spread angle. This suggests a significant 
algorithm dependency on the spatial orientation and movement of the TOA sensors. On 
the whole, a large spread angle between the sensors is preferable. This brings the sensor 
paths closer to orthogonal directions, boosting the potential spatial resolution of the TOA 


measurements. 
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It should be mentioned that an RMSE on the order of 10* m renders the emitter 
estimate effectively useless for any high-accuracy geolocation applications. We suspect 
that the variability in the TOA measurements due to ionospheric propagation is simply 
too large for accurate estimation. For beacon signal propagation over distances of more 
than 19,100 km, small measurement deviations of 9.11 ns are amplified significantly. 
This speculation is validated by the drastic reduction in RMSE for the low-noise trials 
with 10% of the time delay standard deviation. This low-noise condition may be the 
result of atmospheric conditions, such as nighttime operation, or may be the result of time 
delay variation compensation. Systems such as the GNSS use the secondary L2 band at 
1227.6 MHz to boost the accuracy of the TOA measurement [23]. The overall effect is a 
reduction in the random nature of the TOA measurement, corresponding to a lower time 
delay standard deviation and lower RMSE. Such compensation techniques may be 
incorporated into this model in the future using forecast models such as those hosted by 
the CCMC to predict the true time delay characteristics of the ionosphere [18]. Rather 
than aggregating all ionospheric effects simultaneously in a single worst-case random 
variable, the system could be modified to adapt and dynamically calibrate in response to 


current ionospheric conditions. 


An inherent source of error in these results follows from the discretized 
parametric solution strategy used. Each TOA circle was discretized below the set 
maximum altitude h,,,, =5,000 m using 20,000 sample points. As with any numerical 
solution, we have a small-scale variation due to quantization and rounding errors. We 
consider this as a minor effect in the above satellite sensor results in light of a final 


source of uncertainty. 


A final consideration that applies to all results included in this research is that the 
RMSE is approximated here using ten independent trials. Ideally, to approximate the 
expected value operation embedded within the RMSE definition (see Equations (108) and 
(109)), something on the order of 10° or 10° independent trials should be used. Here, 
due to the practical limitations of long compute times, we were only able to simulate ten 


trials. Practically, this implies that the RMSE traces given in Figure 27 contain some 
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undesired variation in the data due to the limited sample size of this stochastic process. 
This is the likely cause of the small-scale variations between data points. If more 
computational resources—including both hardware and time—were available to increase 
the number of iterations as well as decrease the spread angle step size, we anticipate that 
curves above would be smoother and lack the high-frequency deviations seen between 


adjacent data points. 


b. UAV Sensors 


The RMSE for the synthetic aperture TOA algorithm applied to UAV sensors is 


plotted against sensor velocity spread angle @ in Figure 28. 
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Figure 28. | RMSE 1s plotted against sensor velocity angular spread a@ for UAV sensors. 
Dashed traces are linear interpolations between adjacent data points. 
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The simulation results for UAV sensors maintain similar characteristics with the 
satellite sensor results. Both contain uncertainties arising from inherent quantization 
errors and the practical computational limits of the simulation. On the whole, as with the 
satellites, a larger spread angle @ corresponds to a lower RMSE. The one notable 
anomaly is at @=0°, which corresponds to the case when the sensors travel parallel to 
one another. It is unclear what part of the synthetic aperture TOA algorithm would favor 
system performance at such a condition over and against the case with @ =10°. Further 


investigation of this phenomenon is needed. 


Again, we see a pronounced drop in the RMSE for the low-noise condition. For 
the UAV scenario, we observe an RMSE below 50 m for all spread angles. At a 
minimum, the RMSE drops to 12.0 m when @ =90°. Even without specific time delay 
variation compensation measures, geolocation estimates accurate to within 12.0 m are 
very useful. Certainly for applications such as personal navigation or altitude estimation, 


this accuracy rating is sufficient. 


2. Aperture Size 


We next consider the performance of the synthetic aperture algorithm against the 
size of the aperture. The number of aperture samples is varied from N =2—the 
minimum requirement to generate an emitter estimate—up to a maximum aperture size 
N=N,x- AS the number of aperture samples N increases, the time interval between 
samples AZ;, is reduced so as to hold the swath angle traced by the sensors constant. 


Mathematically, we may write 


(113) 


SA,min ? 


iy eae 
AT,, =| —=__|AT, 
SA [As ] 


where AT,, i, 18 the minimum aperture time interval associated with the largest number 


of aperture samples N,,,,. A diagram illustrating the variable synthetic aperture size is 


given in Figure 29. In the diagram, black dots represent the sample times. 
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Figure 29. The synthetic aperture time interval A7,, decreases with increasing synthetic 
aperture size N . Black dots represent sample times. 


As before, 20,000 sample points are used to discretize each TOA solution circle. 
An altitude constraint is applied for the sake of estimate resolution and reduced 
computation time. The low-noise condition again corresponds to a time delay standard 
deviation equal to 10% of the high-noise condition value o,. Other pertinent simulation 


parameters are listed in Table 2. 


Table 2. Aperture size simulation parameters are given for the 
synthetic aperture TOA algorithm. 


Parameter (units) Satellite Sensors UAV Sensors 
19.10 10° 3.048 x 10° 


3.88910" 51.44 
a (degrees) 45 45 


Nx (Samples) 50 50 
N (samples) variable on [2,15] variable on [2,15] 


AT 5 és (8) 64.3 0.778 
1.750x10° 280.1 
1000 100 
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a. Satellite Sensors 


The RMSE for the synthetic aperture TOA algorithm applied to satellite sensors is 
plotted against synthetic aperture size N in Figure 30. Note that we do not plot results 
for synthetic aperture size N <2. Recall that the minimum requirement for the synthetic 
aperture algorithm is to have more than one measurement in time. With only one sample 
point in the aperture, we have a single TOA circle at that given time. We know the 
emitter lies somewhere on the circle in the noiseless case, or near the circle in the noisy 
case, but we cannot refine the search beyond this. There is insufficient data. Thus, we 


begin all plots with the synthetic aperture size with N =2. 
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Figure 30. | RMSE is plotted against synthetic aperture size N for satellite sensors. 
Dashed traces are linear interpolations between adjacent data points. 


As when we plotted the RMSE versus spread angle, we must mention the very 


large amplitude RMSE observed when using satellite sensors. For the sake of brevity, we 
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will not reiterate the discussion on mitigating noise effects using time delay variation 
compensation techniques. Instead, let us highlight the overall trend clearly seen with 
increasing synthetic aperture size. For the high-noise condition, we observe a minimum 
RMSE of 1166 m for N = 50; for the low-noise condition, the minimum RMSE is 125.8 
m for N = 44. On the whole, a larger synthetic aperture corresponds to lower RMSE and 
a better emitter location estimate. Intuitively, this makes sense. If we collect more data 
concerning a random process within a given time interval, the refined aggregate estimate 
should improve. There are some small-scale deviations from this overall trend; however, 
due to the small number of iterations, we do not recommend significant conclusions be 
drawn from such small-scale variations. Significantly more than ten iterations of the 
program are necessary before drawing small-scale inferences from the data. Still, the 


overall trend seems valid. 


b. UAV Sensors 


The RMSE for the synthetic aperture TOA algorithm applied to UAV sensors is 
plotted against synthetic aperture size N in Figure 31. Minimum RMSE values of 55.8 
m (N =16) and 2.89 m (N = 29 ) were observed for the high- and low-noise conditions, 
respectively. The most notable feature of these results is the anomalous trend seen for the 
high-noise condition, namely that the RMSE seems to increase with the synthetic 
aperture size beyond the local minimum at N =16. This is certainly a counterintuitive 
finding, especially since the trend is reversed for the low-noise condition. It appears that 
for low-altitude UAV sensors, increasing the synthetic aperture size allows a more 
significant degradation of the emitter location estimate due to high noise. We speculate 
that there is a point of diminishing return with the synthetic aperture principle. If, for a 
given altitude, the time delay variation is below a certain threshold standard deviation 
value, the synthetic aperture algorithm will indeed improve the estimate (see the low- 
noise data). Alternatively, if the time delay variation exceeds that standard deviation 
threshold, it seems that increasing the number of aperture samples actually degrades the 
accuracy of the emitter location estimate. Granted, these RMSE results are 


approximations based on ten iterations, and a program with a larger sample size will 
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reveal the true theoretical trend for RMSE versus aperture size. Still, we conclude that 


this is a very significant result. 
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Figure 31. RMSE is plotted against synthetic aperture size N for UAV sensors. Dashed 
traces are linear interpolations between adjacent data points. 


D. SUMMARY 


Overall, the synthetic aperture concept was successful in enabling geolocation 
estimates to be formed using only two airborne sensors. For satellite sensors, future 
implementations of this synthetic aperture TOA algorithm require some form of time 
delay variation compensation techniques in order to deliver results with a useful degree 
of accuracy. For UAV sensors, the RMSE performance under low-noise conditions was 
quite acceptable. Overall, the system accuracy tends to improve when the sensors travel 
with a spread angle near @ = 90°, with the curious anomaly of high-performance with 
parallel sensor motion. For normal to low noise conditions, increasing the number of 
aperture samples N results in improved RMSE performance. An interesting exception to 


7S 


this trend arose for high-noise conditions with UAV sensors. Here, it seems that 
increasing the number of aperture samples actually worsened the emitter location 


estimate. 
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VII. SYNTHETIC APERTURE: TDOA-FDOA FUSION 


A. OVERVIEW 


With the TOA algorithm, the three-dimensional intersection curve of two emitter 
surfaces at a given time is always a circle that may be parametrically defined. When we 
move to consider the fusion of one TDOA and one FDOA emitter surface at a given time, 
we are not so fortunate as to have a simply-defined analytical solution. Recall that for 
TOA, we generate one emitter surface per sensor, whereas for TDOA and FDOA, we 
only generate one emitter surface per two sensors. As such, with two airborne sensors 
equipped with TDOA and FDOA systems, we can calculate a maximum of two emitter 
surfaces at any given point in time. The TDOA governing equation (see Equation 2) and 
the FDOA governing equation (see Equation 41) are very difficult to combine. Most of 
this difficulty arrives from the complexity of the three-dimensional FDOA emitter 
surface. With TDOA, the solution is always one half of a hyperboloid of two sheets; this 
is relatively simple to define analytically. Unfortunately, the FDOA surface is defined by 
the solution to a first-order differential equation; the governing FDOA equation relates 
the difference in arriving frequency with the vector forms of each sensor velocity—the 
first-order derivative of the sensor position. The practical consequence is that we require 
a very different solution strategy for the synthetic aperture TDOA-FDOA fusion 
algorithm than was used for the synthetic aperture TOA algorithm. 


B. SOLUTION STRATEGY 


Since an analytical solution is too complicated, an exhaustive search solution is 
selected as an alternative. Whereas the parametric TOA solution considers the discretized 
geolocation solution curve derived from the TOA measurements and calculates the 
spatial coordinates of the estimate, this TDOA-FDOA fusion solution considers a 
discretized region of spatial coordinates and determines which coordinates best satisfy 
the geolocation measurements. This strategy may be seen as iteratively reverse- 


engineering the problem. 


TT 


iP Bounding Box Definition 


To begin, we select a region of three-dimensional space where we expect the 
emitter to be located. For the sake of computational implementation and simplicity, 
consider this region to be a cubic volume centered around the initial seed coordinates 
S (Ags Vesey} with side length equal to 2g. We shall iteratively search this volume to 
find the most-likely emitter location. Since this region is cubic, it shall be referred to as 


the algorithm’s bounding box. A basic diagram of the bounding box is given in Figure 32 


when the point of origin is selected for the seed coordinates. 


bounding box 





Figure 32. A bounding box is centered about the coordinate system origin. 


Zs Discretization Technique 


Next, the bounding box is discretized into a three-dimensional constellation of 
sample points. For this algorithm, we apply a uniform sampling density function, 
meaning we sample the bounding box in each dimension with K points at a constant 


spatial interval equal to (2g)/(K-1). The result is a uniformly-dense constellation 


consisting of K° total sample points. We map the coordinates (% y,Z) of each 
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constellation point to a set of matrix indices (i, j,k), where i, j, and k are bounded on 


the range [1, K] ‘ 


3 TDOA and FDOA Metric 

The algorithm then iterates through each point in the bounding box matrix for 
each time step. For each sample point, we compute the TDOA and FDOA values that 
would theoretically be received by sensors s, and s, if an emitter were located at that 
point. The result is two KxKxkK matrices of TDOA and FDOA measurements. We 
then apply a normalized Gaussian metric to each of these matrices. Based upon the 


TDOA measurement matrix M,, we populate the TDOA metric matrix M , for all 
indices (LiKE 


-(M,.(i, j,k)—t) 


2 (Gmai ; 


(114) 





M, (i, j,k) =exp 


where T,, is the actual received noisy TDOA measurement and 0;,,, is an empirically- 
defined parameter. The metric is thus bounded in the range (0,1] . For a relatively small 
Orpo, Value, the metric value falls off rapidly toward zero as M, (i, j,k) deviates from 
the TDOA measurement T,,. Conversely, a relatively large 0O,,,, value assigns metric 
values close to unity even for substantial deviations of M,(i, j,k) from T,. In this 


algorithm, we select 


Trine = (35 |[mex (ot, )-min(,)) (115) 

The 1/20 factor is a tunable decision parameter. For applications with a fine spatial 

sampling interval, a smaller factor (<1/20) is recommended and vice versa for 
applications with coarse sampling. 

A similar metric matrix M y 18 populated based on the FDOA measurement 


matrix M ,. as a function of indices (i, j,k): 
719 


-(M (i, j.k)-A fy) 





M , (i, j,k) =exp 5 : (116) 
2 (e570) 
where A f,, is the received noisy FDOA measurement and O;,, is chosen to be 
1 . 
Srnos =| 5p |Emx(Mp)~min(M)} (117) 


Each element in the FDOA metric matrix M , 1s bounded on the range (0,1] . 


A composite metric M , 18 defined as the sum of the TDOA and FDOA metrics 


for a given time step 1. Over a total of N synthetic aperture samples, we express the 
aggregate composite metric as 


M=>M,=>(M,,+M,,) (118) 


n=l n=l 
Preserving our double-primed notation from Chapter VI, the estimated emitter position is 
chosen to be the coordinates (x",y",z") that are mapped to the matrix indices (i”, j",k") 
at which the aggregate composite metric is maximized. Mathematically, the estimated 


emitter position vector r= x"x + y"y+z"Z is chosen such that 
M (i, j,k") =max| M (i, j,k) |. (119) 


4. RMS Error 


For the synthetic aperture TDOA-FDOA fusion algorithm, we again use the 
RMSE as our performance metric. The RMSE is defined identically as in Equations 
(106)-(109) in Chapter VI regarding the synthetic aperture TOA algorithm. To avoid 


redundancy, we shall not repeat the same discussion. 


C. RESULTS 

The following results are generated based on the same scenario as considered in 
Chapter VI, with the following modifications. The stationary emitter e generates a 
combined TDOA and FDOA beacon signal that is recetved by the sensors s, and s,, 


which both are equipped with TDOA and FDOA receiver systems. The TDOA 
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measurement is randomly affected according to a time delay standard deviation V20.,: 
Additionally, the deterministic FDOA measurement is affected by additive zero-mean 


Gaussian noise with standard deviation \2o, > as discussed in Chapter V. 


1. Sensor Velocity Angular Spread 


For RMSE simulations against the sensor velocity angular spread a , we vary the 
sensor velocity v, as detailed in our Chapter VI discussion. The primary difference 
between the TOA and TDOA-FDOA algorithms here is that we do not use the maximum 


altitude h,,,, constraint. Instead, the TDOA-FDAO search requires delineation of the 


jax. 


bounding box defined in terms of the seed coordinates S(x,,y,,z;) and the side half- 


length g. For low-noise trials, we use 25% of the high-noise g value listed in Table 3, 
since it is assumed that the low-noise TDOA and FDOA measurements are more 
accurate. We thus limit the search to a smaller region to boost resolution. Pertinent 


simulation parameters are listed in Table 3 for both satellite and UAV sensors. 


Table 3. Angular spread simulation parameters are given for the 
synthetic aperture TDOA-FDOA fusion algorithm. 


Satellite Sensors UAV Sensors 
19.10 10° 3.048 x 10° 


3.889 <10" 51.44 








a (degrees) variable on [0,90] variable on [0,90] 
N (samples) is 5 





AT\, (8) 225 2.72 





1.750x10° 280.1 





1000 100 








1371 0.412 





Ab 
S (x50 Ze) Ga) (0,0,0) (0,0,0) 





4x10° 2x10" 
K (samples) 151 151 
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The values selected for g in these trials were chosen based on the order of 
magnitude of the RMSE observed for similar TOA simulations (see Chapter VI, Section 
C). If a tracking algorithm is implemented along with this system that is able to predict 
future emitter position and if it is known that the emitter is in a smaller region, the 
bounding box specified by half side-length g may be reduced in size. As well, if more 
computational resources are available, the number of one-dimensional bounding box 


samples K may also be increased. 


a. Satellite Sensors 


We have plotted the RMSE of the synthetic aperture TDOA-FDOA fusion 


algorithm for satellite sensors against sensor velocity spread angle @ in Figure 33. 
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Figure 33. | RMSE is plotted against sensor velocity angular spread q@ for satellite 
sensors. Dashed traces are linear interpolations between adjacent data points. 
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Regarding these results, first note that the RMSE for both high-noise and low- 
noise conditions is similar in magnitude to the associated bounding box side half-length 
g. Recall that for low-noise, the bounding box is restricted to 25% of the high-noise 
linear dimension size. This suggests that there is significant variation regarding the 
emitter estimate within the bounding box with limited accuracy. We attribute this poor 


performance to the extremely large quantization errors associated with the bounding box 


search. The spatial sampling interval was set to be (2 g) / (K 1) in each dimension. For 


the satellite sensors that function at extremely high altitudes, the associated bounding 
box, and thus g, must be very large to accommodate the geometry. As g becomes very 
large, the spatial sampling interval likewise increases. Similar to the Nyquist sampling 
criterion in digital communications, if our sampling interval becomes too large, we lose 
the ability to discern high-frequency variations in the associated data. In this case, as the 
spatial sampling interval increases, the discretized TDOA and FDOA measurement 


matrices M, and M ,. lose the ability to discern small-scale variations in the TDOA and 


FDOA equations. The attendant consequence is that spatial resolution decays. As spatial 
resolution decreases, this quantization error likely masks any otherwise prominent trends 


in the data, as is verified by the RMSE measurements presented in Figure 33. 


The obvious caveat to this degradation of spatial resolution is that we may 
increase the number of sample points K to counteract the effects of a large bounding 


box. The practical drawback is that the iterative synthetic aperture algorithm execution 


time is on the order of O(NK : ); we exhaustively compute the associated M, and M 


matrix values for each of K° matrix entries and then iterate the process for each of the 
N synthetic aperture samples. In short, the computation time is prohibitive for high- 
resolution searches of large bounding boxes. Note that for this simulation, K was chosen 
as 151. This may seem curiously small, but consider that the simulation was run within a 
shell script which iterates the algorithm for ten different spread angles, each of which is 
iterated ten times to compute approximate the RMSE value. Altogether, it required over 
seven hours of computation time to generate each of the high-noise and low-noise data 
sets in Figure 33. We see here an immediate problem with a non-analytical solution 
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strategy for large-scale geolocation problems; the computation time rapidly becomes the 
primary limiting factor. Some potential improvements on this point are discussed below 


following the data presentation. 


The significant impact of the quantization error on these results tempers the 
conclusions that might be drawn from the data regarding the impact of sensor velocity 
spread angle. Still, though mitigated in its accuracy, we see a similar trend as with the 
synthetic aperture TOA algorithm. There is slight improvement in RMSE performance as 


the spread angle closes toward a@ =90°. 


b. UAV Sensors 


The UAV sensor results for the similar case of RMSE as a function of spread 


angle is given in Figure 34. 
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Figure 34. | RMSE is plotted against sensor velocity angular spread a@ for UAV sensors. 
Dashed traces are linear interpolations between adjacent data points. 
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With a smaller system geometry, we note that with UAV sensors the RMSE data 
falls further below the g threshold (2,000 m for high-noise and 500 m for low-noise). For 
the low-noise condition, we see a minimum RMSE of 210.8 m at the spread angle of 
a = 20°; this is less than 50% of the bounding box side half-length and suggests that the 
algorithm is indeed gaining some information on the location of the emitter within the 
bounding box. This is an encouraging result, yet the algorithm is still seriously limited by 
the spatial resolution of the bounding box. Unfortunately, the data does not show 
significant trends relating to sensor velocity spread angle a . We hesitate to draw further 


conclusions based on small-scale variations due to our sampling error. 


2. Aperture Size 
Here again, we assess the algorithm performance against the synthetic aperture 
size N as detailed in Chapter VI. We vary the synthetic aperture time interval A7,, in 


accordance with N to maintain a constant swath angle. All pertinent simulation 


parameters used to generate the following results are presented in Table 4. 


Table 4. Aperture size simulation parameters are given for the 
synthetic aperture TDOA-FDOA fusion algorithm. 


Satellite Sensors UAV Sensors 
19.10 10° 3.048 x 10° 








3.889 <10" 51.44 





a (degrees) 45 45 





N (samples) variable on [2,15] variable on [2,15] 








1,750x10° 280.1 





1000 100 











N.ux (samples) 15 15 
AT 4 min (S) 64.3 0.778 
S( zs) 


(000 (000 





4x10' 2x10" 











K (samples) 151 151 
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a. Satellite Sensors 


As with the satellite sensor RMSE results as a function of spread angle, the data 
for RMSE versus synthetic aperture size presented in Figure 35 is dominated by 
quantization error due to sampling and computational resource limitations. For a point of 
reference, the total execution time for each of the high- and low-noise data sets shown in 
Figure 35 was over 16 hours. Accordingly, we do not see the appreciable improvement in 
estimate accuracy that was achieved for the analogous case with the synthetic aperture 
TOA algorithm (see Figure 30). Still, there are interesting features in the data that merit 


discussion. 
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Figure 35. | RMSE is plotted against synthetic aperture size N for satellite sensors. 
Dashed traces are linear interpolations between adjacent data points. 


Note the nearly point-for-point match in trend between the low-noise and high- 


noise data sets. We observe local maxima at N =2,6,9, and 14 and local minima at 
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N =8,10,13, and 15. Intuitively, we expect the RMSE to monotonically decrease for 


larger N ; as we include more data in the synthetic aperture, the estimate accuracy should 
improve. Since all data was collected based on a series of independent trials, we cannot 
easily attribute this repeated oscillatory behavior to stochastic causes. Rather, we 
speculate that this predictability is a vestige of quantization errors associated with the 
discrete search algorithm. Since the emitter estimate is chosen from a finite set of spatial 
coordinates within the bounding box, some of the small-scale random nature of the noisy 
TDOA and FDOA measurements is absorbed by the rounding process associated with 


quantization. 


b. UAV Sensors 


The RMSE results for UAV sensors are plotted against the aperture size in 
Figure 36. 
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Figure 36. RMSE is plotted against synthetic aperture size N for UAV sensors. Dashed 
traces are linear interpolations between adjacent data points. 
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Here, as before, we see some improvement in the RMSE as compared with the 
bounding box side half-length g; however, the overall absence of a discernable trend 
suggests minimal correlation between the RMSE and aperture size for the algorithm as 
described in this chapter. We had expected to see some improvement with larger aperture 
size, as was seen for the low-noise UAV sensor scenario with the synthetic aperture TOA 
algorithm. Unfortunately, the sampling limitations likely masked any such trends. The 
primary implication to draw from these TDOA-FDOA fusion results is that the sampling 
constraints must be addressed if the algorithm performance is to be refined to a practical 


level. 


D. POTENTIAL IMPROVEMENTS 


To that end, we offer some ideas for potential improvements of the TDOA-FDOA 
fusion algorithm that we could not implement for want of time. One technique would be 


to manipulate the quantization method used on the bounding box to group more samples 


near the seed coordinates S (ates Vioks) and fewer near the perimeter of the box. If the 


seed value is near the emitter, this will boost the spatial resolution of the iterative search 
in the local region surrounding the seed. This effectively focuses the samples within a 
subregion of the bounding box. The drawback to this technique is if the seed coordinate is 
far from the true emitter position, the sample density surrounding the emitter may 
actually become sparser than with a uniform sample density. Such a sample density is 
given as a red trace in Figure 37. For comparison, the uniform sample density used in the 


above model is given as a black trace. 


A more sophisticated alternative is to use an adaptive sampling technique that 
dynamically redefines the sampling density function at each time step in the aperture. 
Based on received measurements, the algorithm might focus the sample density on 
different regions of the bounding box. Such a technique may improve the spatial 
resolution and mitigate quantization errors, but it alone will not reduce the computational 


demands of the iterative search algorithm. 
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Figure 37. A manipulated sample density function (red trace) condenses more samples 
near the seed coordinates than a uniform density function (black trace). 


A second possible improvement is to implement an adaptive telescopic zoom that 
progressively tightens the bounding box to a smaller and smaller region of interest. With 
each telescopic iteration, the spatial resolution of the algorithm improves. A two- 
dimensional illustration of such an algorithm is given in Figure 38. For each iteration in 
this example, the algorithm selects a subregion with 1/9 the total area of the current 


bounding box. Still, this does not ease the demand on computational resources. 





Figure 38. An adaptive zoom technique iteratively reduces the size of the bounding box. 


To improve the algorithm in regards to execution time, an effort must be made to 


apply loop-flattening techniques. A program with execution time on the order of 
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O(NK *) is extremely difficult to scale. If some inherent matrix operations native to the 


MATLAB” software library that was used in this research can be applied to flatten one of 


the internal loops of the iterative search, the runtime might be reduced to the order of 


O(NK als This would likely yield significant reduction in runtime, allowing smaller 


sampling intervals and better spatial resolution. In broad terms, the algorithm’s 
computational efficiency must be improved, regardless of other modifications, if it is to 


be implemented as a viable real-time geolocation algorithm. 


E. SUMMARY 


Unlike the synthetic aperture TOA algorithm, we did not consider an analytical 
solution to the synthetic aperture TDOA-FDOA fusion problem. Rather, we adopted an 
exhaustive search solution strategy within a certain bounding box region. Results from 
our simulations suggested that the quantization errors associated with our sampling 
constraints were too large to yield useful results for high-accuracy geolocation. As a 
result, we offered several options for future improvements to the algorithm to help 
overcome the limitations due to excessive execution times. Altogether, the complicated 
nature of the TDOA and FDOA geolocation equations necessitated a computationally- 
intensive solution strategy which exhausted our computational resources before the 


algorithm could converge to a useful solution. 
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VU. CONCLUSION 


A. RESEARCH CONTRIBUTIONS 


The initial goal of this research was to investigate the feasibility of two-sensor 
synthetic aperture geolocation techniques. The driving motivation was to reduce the 
physical resource requirements for successful geolocation from three or four airborne 
sensors down to only two. We speculated that the use of a synthetic aperture principle 
similar to that used in radar, sonar, and medical imaging, might be helpful. To begin, we 
laid the groundwork of establishing the theoretical two- and three-dimensional solutions 
for the TOA, TDOA, and FDOA geolocation methods. When possible, analytical 
solutions were given to define the emitter curves and surfaces resulting from a given 


geolocation estimate. 


Subsequently, three separate noise models were developed for TOA, TDOA, and 
FDOA measurements for two distinct cases. First, time delay and Doppler shift variation 
parameters were determined based on the assumption that ionospheric signal degradation 
was the primary mechanism of TOA, TDOA, and FDOA measurement variation about 
the deterministic value calculated based on the geometry of the problem. Ionospheric 
forecast models and measurements were considered in the determination of the noise 
parameters. Second, time delay and Doppler shift variation parameters were estimated for 
a geolocation system onboard a UAV. We turned to literature resources to estimate the 
time delay standard deviation and assumed that the Doppler shift variation in the lower 
atmosphere would mainly result from instrumental uncertainties rather than the 
nonlinearity or inhomogeneity of the propagation medium. Together, these parameters 


contributed to a basic noise model for incorporation in our later simulations. 


Regarding the heart of our work—the synthetic aperture algorithms—we first 
investigated the homogenous algorithm consisting of two TOA measurements per time 
step in the synthetic aperture. Using a parametric solution strategy to find the intersection 
points of subsequent TOA emitter circles, we found that the algorithm requires some 


form of time delay variation compensation to function with reasonable accuracy for 
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satellite sensors. Conversely, for UAV sensors the algorithm has promising levels of 
accuracy, especially in low-noise environments. A primary finding is that the RMSE of 
the emitter estimate tended to decrease with increasing synthetic aperture size, except 
perhaps for particularly high noise levels relative to the system geometry. This helped 
validate the general principle of using a synthetic aperture in time to refine the final 


emitter location estimate. 


Applying the synthetic aperture concept to a heterogeneous combination of 
TDOA and FDOA measurements was far more challenging from a solution strategy 
perspective. We were forced to consider an exhaustive search solution strategy in lieu of 
a straightforward analytical solution as for TOA. The results collected from our 
simulations highlighted the practical limitations of an exhaustive search in three- 
dimensions; the computation time required to search a large spatial region with fine 
resolution can rapidly become unmanageable. We are thus faced with a trade-off between 
computation time and spatial resolution of the algorithm. More focused research is 
required in this area before the synthetic aperture TDOA-FDOA algorithm can be useful 


in real-time applications. 


On the whole, we recommend use of the synthetic aperture TOA algorithm for use 
with small-scale UAV geolocation networks. The mathematical solution for the 
homogeneous TOA problem is a distinct advantage over the heterogeneous TDOA- 
FDOA problem. Yet, if the issues related to sampling intervals in the TDOA-FDOA 
algorithm can be solved or mitigated, then it may be found to have advantages over the 
TOA algorithm with respect to RMSE performance as a function of sensor spatial 


orientation or other system variables. 


B. OPPORTUNITIES FOR FUTURE WORK 


In broad terms, the foreseeable opportunities for future work on this two-sensor 
synthetic aperture geolocation problem may be grouped into three categories. First, a 
more robust consideration of noise effects and associated noise mitigation techniques will 
be helpful. Second, the proposed algorithms offered in this work are heuristic in nature; 


work should be directed toward developing optimal estimators for the TOA and TDOA- 
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FDOA cases. Third, more concentrated research should be directed toward determining 
optimal synthetic aperture configuration. Other opportunities for further work do exist, 


but we have tried to highlight those most prominent areas. 


iF Noise Effects and Mitigation 


Although reasonable for initial consideration, the treatment given to noise 
modeling in Chapter V was simplistic. It is unlikely that all atmospheric effects on a 
geolocation signal can be folded into one single random variable with a high degree of 
fidelity. Modifying the time delay noise models to include more generalized PDFs such 
as the Rayleigh, Rician, and lognormal distributions may offer a more robust time- 


domain model. Experimental data should be collected to help validate such models. 


Certainly, work should be done to incorporate the available ionospheric forecast 
models into geolocation systems. Based on the latitude, longitude, and time of day, we 
can much more accurately predict the time delay and Doppler shift associated with 
ionospheric propagation for satellite-based systems. Moreover, we might consider out-of- 
band signaling with known ground stations to sample the channel. Depending on the real- 
time channel condition, appropriate compensations may be incorporated into the emitter 
estimation algorithm. All this is in an effort to reduce the degree of randomness 


associated with the propagation of an emitter’s beacon signal. 


2. Optimal Estimators 


As previously mentioned, the heuristic algorithms presented herein are not 
optimal. That is, they have not been compared against the CRLB. For resources 
concerning optimal estimators, see [35] and [36]. Optimization research directed toward 
finding an efficient estimator for both the synthetic aperture TOA and TDOA-FDOA 
problems may give further insight into methods to refine the computational efficiency of 
these algorithms. In practical terms, the crux of the issue is being able to implement a 
synthetic aperture algorithm in a real-time system that boasts high spatial resolution and 


manageable computation time. 
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3. Optimal Synthetic Aperture Configuration 


Moving forward with the synthetic aperture geolocation concept, we must 
determine what sensor orientation, velocity angle, swath angle, and aperture sample size 
achieve optimal results for both the TOA and TDOA-FDOA algorithms. Our results 
showed notable variation in system performance with each of these variables; the next 
step is to find an optimal set of all variables to achieve the most consistent, accurate 
systems performance. This information will guide the implementation tactics of any 
future synthetic aperture geolocation system and, as such, is a significant next step in the 


development process. 
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APPENDIX. MISCELLANEOUS DERIVATIONS 


A. THREE-DIMENSIONAL TDOA EQUATION 


The following derivation follows very closely that given in Chapter III for the 
two-dimensional TDOA solution, which is in turn based on the proof given in [16]. Let 


us begin by restating the governing TDOA equation: 














r,-T,|- r, —¥,| = C7): (120) 


Expressing this vector equation in terms of its Cartesian components, we have 
2 2 2 2 2 2 
(72) +(y-y,) +(Z=%) ca (4%) +(y-y,) a Z=Z,) =C oT, (121) 
where the two sensors 5, and S$, are located at the generic positions (x, y,,z,) and 


(x,,,Z,)- For the specific case when sensors $, and $, are positioned at coordinates 


(—d,0,0) and (d,0,0), we have 


(x-d) ty’ +2 —-J(xtd) ty? +2 =e,7y. (122) 
Rearranging and squaring gives 


Zz 


2 
| (say +y'+2'] = (cota + (r+ay+y'+2'] (123) 
Expanding both sides fully gives 


MoS Dae Ry Hee 


= (C5) + 2cyt ME DINE EE ey Pet +x°42dx+d?+y’? +z’. 


We cancel terms and rearrange to get 


Adx—(e)T) =2(cgty, yx" + 2dxtd? + y?+2’. (125) 


In terms of the parameter a=c,r,,/2, (125) simplifies to 


—dx—a’ =ax’ +2dx+d*+y' +z’. (126) 


(124) 





We square to get 
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Gi 42daxta' =a’ (x t2dx+d?+y'4 z); (127) 


Collecting terms and factoring leaves 


(a a \x ay az =a" (d?-a’). (128) 





Dividing through by a’ (d : -a’) leaves 


2 2 2 
x y Z 
a d-a@ d-a 
From here, the equation may be cast in terms of the traditional hyperbolic parameter 


b=~d’—a’ to give 





=1, (129) 


a (130) 


or we may substitute out the intermediate parameter a to give the final three-dimensional 
hyperboloidal TDOA equation 
4x? a 4y? = 42? 
2 2 2 
(cyt) 4d? -(cyty) — 4d* (cyt) 





=1, (131) 


B. PARAMETERIZING A CIRCLE: ORTHOGONAL UNIT VECTORS 


Two orthogonal unit vectors & and F,,are known to satisfy the dot-product 
criterion 


AK 


i-F, =(d,8+0,9+02)-(A,,2+%,,9+h)2)=0. (132) 


This implies that unit vector w lies entirely within the x-y plane. Evaluating the dot 


product gives 


Uu, 21x TU, 2ly =0. (133) 
Rearranging and squaring gives 

a2. 2 PD Ds 

Uu, Diy =u, iy _ (134) 


By virtue of u@ being a unit vector, we know 


=17 +4. =1, (135) 


x y 
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u 














and thus, 
a, =1-i,’. (136) 
Substituting (136) into (134) gives 
(1-a,7)4,7 =4,74,,”. (137) 
We then distribute terms and factor to get 
hic =o Foy (138) 
Retaining the root ambiguity, we have 
~ 2 
fi, =+/—2—.. (139) 
: Diy thy 


Equation 133 may be solved solve for u, such that 








a A ~ 2 a) 
r, r. r r. 

A 2ly |a pied 21y i) 2ly 

i--(# ye le ee ye (140) 
My Diy) Nix thy Diy thay 


Finally, we may express the unit vector @ in the three-dimensional Cartesian 


coordinate system as 


pe p, 7 
els |r a |p (141) 


Nyy 21ly thy 


C. SCALED APERTURE INTERVAL FOR CONSTANT SWATH ANGLE 


Consider the generalized flat-earth diagram for a single geolocation sensor 5s, at 
altitude h passing over an emitter e given in Figure 39. The sensor travels at constant 
speed v, parallel to the ground. Over the time period Tj, associated with the synthetic 


aperture, the sensor traces out an associated swath angle y . 
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~ 





Figure 39. — A geolocation sensor traces out a swath angle y relative to an emitter. 


The length of the aperture traced by the sensor can be written as 
W = 1, Toy. (142) 
If the synthetic aperture is composed of N individual, regularly-spaced measurements, 
then we may write 
Ty =(N -1)AT,y. (143) 
Substituting (143) into (142) gives 
W =v,(N-1)AT,,. (144) 
From basic trigonometric relationships, we see that 


tn( We Cu EL 
2) 2h 2h 





(145) 


We aim to scale the synthetic aperture time interval A7,, such that the swath 


angle is held constant for an arbitrary sensor velocity and altitude. If the swath angle is 
fixed, (145) is likewise a fixed constant. Accordingly, all scaled values (primed notation) 


are valid so long as we satisfy 
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| 


v,(N-1)AT,, 


h 


Jl 


vi (N'—1)AT,, 
h’ 


= constant. 


(146) 


To specifically scale our aperture time interval for satellite sensors to the appropriate 


value for UAV sensors, we simply solve for the scaled synthetic aperture time interval: 


AT,, 


h' 


(¢ 


v,(N -1) 
(N’=1) 
3889 m/s)(5-1) 


h 


Ss 
Uy 
Ss 


r( 





=2.72s8. 


Jar 


| (51.44 m/s)(5—1) 


SA 


{ 
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3,048 m 


19.10x10° m 


(147) 


)(22s s) 
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